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ABSTRACT 

Wc study the behavior of magnetorotational turbulence in shearing box simulations with a radial 
and azimuthal extent up to ten scale heights. Maxwell and Reynolds stresses are found to increase by 
more than a factor two when increasing the box size beyond two scale heights in the radial direction. 
Further increase of the box size has little or no effect on the statistical properties of the turbulence. An 
inverse cascade excites magnetic field structures at the largest scales of the box. The corresponding 
10% variation in the Maxwell stress launches a zonal flow of alternating sub- and super-Keplerian 
velocity. This in turn generates a banded density structure in geostrophic balance between pressure 
and Coriolis forces. We present a simplified model for the appearance of zonal flows, in which stochastic 
forcing by the magnetic tension on short time-scales creates zonal flow structures with life-times of 
several tens of orbits. We experiment with various improved shearing box algorithms to reduce the 
numerical diffusivity introduced by the supersonic shear flow. While a standard finite difference 
advection scheme shows signs of a suppression of turbulent activity near the edges of the box. this 
problem is eliminated by a new method where the Keplerian shear advection is advanced in time by 
interpolation in Fourier space. 

Subject headings: diffusion — hydrodynamics — instabilities — planetary systems: protoplanetary 
disks — solar system: formation — turbulence 



1. INTRODUCTION 

Turbulence forms a ma in pillar of modern acc retion 
disk theory (see review bv lBalbus & Hawlevl[T998h . The 
most promising turbulen ce candidate is the magn etorota- 
tional instabihty (MRI, iBalbus fc Hawlevl[l99il ). which 
renders Keplerian rotation profiles linearly unstable in 
the presence of a magnetic field of moderate strength. 
The non-linear evolution of the MRI is often modeled in a 
local, corotating box representing a small section of a Ke- 
pleria n disk (e.g. iHawlev et al.lll995t iBrandenburg et all 
|1995[) . This allows high resolution studies with periodic 
boundary conditions. However, most shearing box sim- 
ulations have been limited to a narrow radial extent of 
approximately one scale height. 

In this paper we explore shearing box simulations 
with a radial domain up to ten scale heights. In- 
creasing the radial extent helps bridge the gap between 
smaller local sim ulations and global simulations of ac- 
cretion disks fe.g.lArmitagdll998HArlt fc Riidigal l2001l: 
iFromang fc Nelsonll2006t iLvra et al.ll2008n . Since rotat- 
ing shear flo ws are prone to an inverse cascade o f mag- 
netic energy (jBrandenburg fc Subramanianll2005f ). large 
simulation domains are needed to fully capture the sat- 
urated state of the MRI. 
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We find that the statistical properties of magnetorota- 
tional turbulence converge relatively well as we increase 
the radial extent of the simulation box beyond a couple 
of scale heights. However, large shearing boxes do not 
behave simply like multiple copies of moderately sized 
boxes. Axisymmetric pressure bumps, in geostrophic 
balance with a sub-Keplerian/super-Keplcrian zonal flow 
envelope, grow to fill the radial width of the box. A main 
purpose of this paper is to describe the excitation and 
dynamics of such zonal flows. 

The emergence of zonal flows in protoplanetary disks 
can have a profound influence on the growth of terres- 
trial planets and giant planet cores. Planet formation 
beg ins as dust grains collide and grow to ever larger bod- 
ies (ISafronovl [19691 : iDominik et all [20071 : iBlum fc WuTml 
I2008D . Macroscopic solids drift rapidly through the disk, 
due to the head wind from the slightly sub-Keplerian 
gas ([Wcidcnschilling 1977), imposing a severe time-scale 
constraint on the growth from approximately cm-sized 
pebbles to bodies several kilometers in size. A weak vari- 
ation of the rotation profile caused by zonal flows may 
nevertheless suffice to locally stop th e radial dri ft flow 
jWhipplc 1972: Klahr fc Linll200l[: [Haghig hipotu- fc"Bosi 
[2003t IFromang fc Nelsonll2005t iKatoeraLl I2OO80 . Such 
a convergence zone can quickly accumulate rocks and 
boulders, enhancing the local density enough to trigger 
gravitational instabilities which fo rm planetesimals or 
dwarf planets in a very short time (jYoudin fc Shul[200^ : 
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iJohansen et al.ll2006ll2007f) . The reduced radial drift due 
to persistent zonal flows may also facilitate dust growth 
by a longer tim e-scale coagulation-fragmentatio n pro- 
cess (IWeidenschillin g 1997; DuUcmond & Dominik] l2005l : 
iBrauer et al.ll2008al bi: iJohanseu'etini^OSIl . 

Turbulent density fluctuations exert significant grav- 
itational torques on planetesimals and planetary em- 
bryos. This adds a stochastic element to the normally 
inward type-I planetary migration from gravitational in- 
teractions with the gas disk. Migration can be slowed, 
sped u p, or even reversed, especially for lower mass 
planets llNelson fc Papaloizou '2004': Laughlin et al.ll2004l : 
iJohnson et all 120061 : lOishi et al.. .206W . Gravitational 
scattering can also pump planetesimal eccentricities to 
the level where collisions become dest ructive for sma ll 
planetesimals with low surface gravity (jlda et al.l 120081 ). 
The dominant axisymmetric component of the density 
fluctuations that we analyze in this paper exerts no grav- 
itational torques. Thus non-axisymmetric modes ~ per- 
haps in global simulations - must be considered in a fu- 
ture publication for applications to planet migration and 
scattering. 

Zonal flows arise in a diverse range of physical settings. 
The launching mechanism in the examples below differ, 
both from each other and from the problem at hand. A 
common feature is non-linear mode coupling and varia- 
tion in turbulent transport coefficients over large length 
scales. 

Giant planet atmospheres. The banded structure in 
Jupiter's cloud layer is perhaps the most famous mani- 
festation of zonal flows. iBussd ()1976[ ) proposed that these 
flows arise from an inverse cascade of thermally driven 
convection in giant planet interiors. Simulati ons of con- 
vective turbul ence in spheric al shells (e.g., ISun et aTl 
[1993: Hcimpel fc Aurnoull2007l ) reproduce the zonal flows 
in gas and ice giants both qualitatively and quantita- 
tively. 

Torsional oscillations. On top of the differential ro- 
tation profile of the Sun there is a zonal flow (with 3 
m/s amplitude) that mi grates from high to l ow la titudes 
during the solar cycle (jHoward fc Labontd 119801 ). The 
zonal fiow is attributed to the toroidal component of the 
Lorentz force exerted by magnetic fields arisin g in the 
solar dynamo (jSchussleii 1 1 98 ll: I Yoshimur al 1 1 98 It) . 

Laboratory plasmas. Laboratory plasmas often exhibit 
strong zonal flows. These self-organized structures arise 
from a non-local invers e cascade mediated by drift waves. 
The excellent review by [Diamond et all (j2005} ) points out 
the intimate mathematical equivalence between Rossby 
wave dynamics in planetary atmospheres and drift wave 
dynamics in plasmas. 

In magnetorotational accretion disk turbulence we flnd 
that the zonal flow is excited by a large scale variation 
in the Maxwell stress. The differential momentum trans- 
port by the associated magnetic tension launches the sub- 
Keplerian/super-Keplerian zonal flows, which in turn are 
slightly compressive and form pressure bumps with life- 
times approximately equal to the turbulent mixing time- 
scale. The variation in the Maxwell stress is a conse- 
quence of an inverse cascade of magnetic energy to the 
largest scales of the box. 

Practical numerical issues arise in simulations of large 
shearing boxes, since the supersonic orbital shear intro- 
duces significant numerical diffusivity. In simulations 



of magnetorotational turbulence, IJohnson et al.l ()2008f ) 
found a spurious density depression in the radial center 
of the grid. Here the velocity of the shear fiow vanishes, 
minimizing numerical dissipation, and the local increase 
in turbulent stresses leads to a mass transport away from 
the box center. We also find indications of this unphysi- 
cal behavior in our biggest simulations, if we use the stan- 
dard shearing box solver of the Pencil Code. However we 
can essentially eliminate this undesired effect with two 
independent techniques described in the appendices: (1) 
separately solving the known advection by orbital shear 
with high order interpolation or (2) displacing the entire 
grid systematically in the radial direction so that pre- 
ferred locations in the box vary too rapidly to generate 
artificial structures. This gives us confidence that the 
zonal flows we see are not a numerical artifact. 

The structure of the paper is as follows. In ^ we de- 
scribe the shearing box coordinate frame, dissipation op- 
erators and simulation parameters. The properties of 
turbulence in large shearing boxes are presented in ^ 
The following section, 21 is dedicated to describing the 
large scale zonal flows that are ubiquitous to the sim- 
ulations. In fjni we present a simplified model for the 
stochastic excitation of zonal flows. The radial variation 
in Maxwell stress and the inverse cascade of magnetic 
energy is discussed in detail in ^jSl In 33 we validate 
our results by presenting convergence tests, variations of 
the main simulations using different shearing box algo- 
rithms, and the dependence of the zonal flows on the 
order and magnitude of the explicit dissipation. In SjSjwe 
summarize our results briefly and discuss caveats and im- 
plications. The appendices contain a test of the shearing 
box solver (Appendix [X| , a description of a newly de- 
veloped technique to interpolate Keplerian advection in 
wavenumber space (Appendix [B| , and additional tech- 
niques to reduce space-dependent numerical diffusivity 
(Appendix [C]). 

2. MODEL SETUP 

We simulate a local, corotating patch of a Keplerian 
disk in the shearing box approximation. The axes are 
oriented such that x points outwards along the cylindri- 
cal radius, y points along the main rotation direction, 
and z points vertically out of the disk parallel to the 
Keplerian rotation vector fl. The criterion for the valid- 
ity of the shearing box approximation is strictly that the 
box size L should be much smaller than the orbital radius 
r, so that curvature terms may be ignored. In practice 
this implies that the box size measured in scale heights, 
L/ H, must be much smaller than the inverse aspect ra- 
tio of the disk, (L/H) ^ l/{H/r). For typical values 
(H/r) = 0.01 ... 0.1 we get (L/H) < 10 . . . 100. We shaU 
consider boxes with a radial extent of up to ten times the 
scale height, which is still fairly well represented in the 
shearing box approximation. We have chosen the shear- 
ing box approximation over global disk simulations be- 
cause of the absence of boundaries and because it makes 
it possible to study convergence in resolution and box 
size from a relatively well understood coordinate frame. 

The ideal MHD equation of motion, for the velocity 
field u relative to the linearized Keplerian shear, is 
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+ ^JxB--VP + f,{u,p). (1) 
P P 

The left hand side of the equation contains the advection 
by the perturbed velocity u and by the Keplerian shear 

flow Uy^By, where Uj,"' ~ — (3/2)!?x. The right hand side 
contains the usual terms for Coriolis force, vertical com- 
ponent of the gravity of the central object, Lorentz force, 
pressure gradient force and an explicit viscosity term 
constructed to dissipate the kinetic energy released by 
Reynolds and Maxwell stresses (see ^2.1l for details). The 
magnetic field is calculated from the magnetic vector po- 
tential through S = V X A, and the current density is 
calculated through Ampere's law J = fj.^^'V x (V x A). 

The magnetic vector potential A is evolved through 
the uncurled induction equation 



dA 

'dt 



,(0) 



dA 

dy 



u X B 



f,iA). (2) 



The terms on the right hand side are the electromotive 
force and an explicit magnetic stretching term represent- 
ing the creation of from Ay (conversely creation of 
By from Bx) by the Keplerian shear. Explicit resistiv- 
ity is taken into account through the term (see W2.1\\ . 
For simulations with a weak imposed vertical field, we 
add a constant field {Bz)ez to the magnetic field B ap- 
pearing in the electromotive force in equation ^ and 
in the Lorentz force in equation ([T]). In the analysis of 
our results we will seek to understand the evolution of the 
magnetic field via the action of individual terms from the 
induction equation for so we give it here for reference. 



{B-V)u--QBxey- 



BVu. 

(3) 

The advection terms appear in their usual form on the 
left hand side. The terms on the right hand side of equa- 
tion ([3]) are the internal stretching term, the Keplerian 
stretching term, and the compression term. 

Finally the mass density p is evolved through the con- 
tinuity equation 



dp ^ 



y dy 



ui:'> ^ = -pV -u + fu {p) ■ (4) 



The last term on the right hand side is an explicit mass 
diffusion term (presented in ^2.ip . We use an isothermal 
equation of state, P = c^p, to connect the pressure P 
and the density p. The soundspeed Cg is assumed to be 
constant. 

As a numerical solver we use the 6th order symmet- 
ric finite difference code Pencil Code^ . In Appendix [X] 
we test the standard shearing box algorithm of the Pen- 
cil Code against a time -dependent analytical solution 
(jBalbus fc HawlevI 12006( 1 for a shearing wave and find 
excellent agreement. We also compare the evolution of a 
low amplitude, magnetized shear wave to a numerical in- 
tegration of the linearized evolution equations. However, 
while the tests in Appendix [A] give confidence that the 
Pencil Code solves the shearing box equations correctly, 

^ The code is accessible for download at 

http : //www.nordlta. org/sof tware/pencil-code/ | 

see Brandenburg (2U()3j tor dlSafls onT;he numerical algorithm of 

the Pencil Code. 



they are not sufficient to fully validate our non-linear sim- 
ulations. For this we will present the main results using 
several independent shearing box algorithms. Normally 
we use the same finite difference advection scheme for 
the Keplerian shear advection as for any other advection 
term, but we have also implemented the Keplerian advec- 
tion as an interpolated shift of the dynamical variables 
(simi lar to iMassetl 120001: iGammid 1200 It I Johnson et"an 
|2008[ ). This has the advantage that the space-dependence 
of the numerical diffusivity is removed and that the time- 
step constraint connected with the Keplerian shear flow 
is eased. Our new Shear Advection by Fourier Interpola- 
tion (SAFI) scheme, described in detail in Appendix [Bl 
differs from other shear advection schemes in that inter- 
polation is done in Fourier space, yielding essentially per- 
fect interpolation of sufficiently smooth flows when shift- 
ing a mode by any arbitrary fraction of a grid cell. Two 
additional algorithms for reducing the space-dependence 
of the numerical dissipation are presented in Appendix [Cl 

2.1. Dissipation 

The symmetric finite difference scheme of the Pencil 
Code has very little intrinsic amplitu de error (i.e. numer i- 
cal diffusion) in the advection terms ()Brandenburg|2003f) . 
Therefore we must add explicit dissipation terms to the 
evolution equations in order to get rid of the kinetic and 
magnetic energy, released from the gravitational poten- 
tial by Reynolds and Maxwell stresses, on the smallest 
scales of the grid. 

2.1.1. Viscosity 
The viscosity term f^^ appears in full generality as 



u + 2{S^^^ • Vlnp) 
V^M+ (S'f^' • Vlnp) 



-j/sh [VV • u 



(V-w)(Vlnp)] + (Vi.,h)V-M. (5) 



The viscosity includes regular Navier-Stokes viscosity 
with constant coefficient vi, hyperviscosity with con- 
stant coefficient and shock viscosity with a variable 
coefficient Vgh- In general we include only a subset of 
the three types of viscosity, representing both direct nu- 
merical simulations (DNS) and combined hyperviscos- 
ity/shock viscosity solution of the equation of motion. 
In the following paragraphs we discuss the three types of 
viscosity in more detail. 

Navier-Stokes viscosity (i^i): Since the density is not 
nearly constant in the stratified models, we must add a 
density dependent term 2(5'^^-' • Vlnp) to equation ^ 
in order to conserve momentum. Here the traceless ratc- 
of-strain tensor S''^-' is 



o(l) _ 

- 2 



1 



dxj dxi 3 



(6) 



Hyperviscosity {v^): A simplified third order rate-of 
strain tensor S^^"^ is here defined as 



:(3) 



dx^j 



(7) 



The high order Laplacian in equation ([5]) is expanded 
as V*^ = d'^/dx^ + d^/dy^ + d^/dz^. This form of the 
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TABLE 1 

Run Parameters 



Run 




X Ny X 






m 




Order 


Strat. 




Shear 


At 


S 


1.32 X 1.32 X 5.28 


32 X 32 X 128 


6.0 X 10" 


-10 


6.0 X IQ- 


- 10 


3 


Yes 


0.0 


FDA 


100 


iVl 


Z.D4 X Z.D4 X i3.2o 


04 X 04 X izo 


D.O X lU 


-10 


n \y 1 n— 
O.U X iU 


-10 


Q 
O 


Yes 


n n 
U.U 


r UA 


lUU 


L 


5.28 X 5.28 X 5.28 


128 X 128 X 128 


6.0 X 10" 


-10 


6.0 X 10" 


-10 


3 


Yes 


0.0 


SAFI 


100 


H 


10.56 X 10. 5d X 5.28 


256 X 256 X 128 


6.0 X 10~ 


-10 


6.0 X 10" 


-10 


3 


Yes 


0.0 


SAFI 


100 


M_r2 


2.64 X 2.d4 X 5.28 


128 X 128 X 256 


2.0 X 10~ 


-11 


2.0 X 10" 


-11 


3 


Yes 


0.0 


FDA 


100 


L_nogz 


5.28 X 5.28 X 5.28 


128 X 128 X 128 


6.0 X 10- 


-10 


6.0 X 10" 


-10 


3 


No 


0.0 


SAFI 


100 


M_Bz_r0.5 


2.64 X 2.64 X 5.28 


32 X 32 X 64 


2.0 X 10" 


-8 


2.0 X 10- 


-8 


3 


Yes 


0.01 


FDA 


100 


M.Bz 


2.64 X 2.64 X 5.28 


64 X 64 X 128 


6.0 X 10" 


-10 


6.0 X 10" 


-10 


3 


Yes 


0.01 


FDA 


100 


M_Bz_r2 


2.64 X 2.64 X 5.28 


128 X 128 X 256 


2.0 X 10" 


-11 


2.0 X 10" 


-11 


3 


Yes 


0.01 


FDA 


50 


MS 


2.64 X 2.64 X 1.32 


64 X 64 X 32 


6.0 X 10- 


-10 


6.0 X 10" 


-10 


3 


No 


0.0 


SRD 


100 


MSj-2 


2.64 X 2.64 X 1.32 


128 X 128 X 64 


2.0 X 10- 


-11 


2.0 X 10- 


-11 


3 


No 


0.0 


FDA 


100 


MSj-2_fix 


2.64 X 2.64 X 1.32 


128 X 128 X 64 


6.0 X 10- 


-10 


6.0 X 10" 


-10 


3 


No 


0.0 


SRD 


100 


SS.r4_nul 


1.32 X 1.32 X 1.32 


128 X 128 X 128 


3.0 X 10- 


-4 


8.0 X 10- 


-5 


1 


No 


0.0 


FDA 


40 


SS.r8.nul 


1.32 X 1.32 X 1.32 


256 X 256 X 256 


3.0 X 10" 


-4 


8.0 X 10" 


-5 


1 


No 


0.0 


FDA 


25 


MS_r4_nul 


2.64 X 2.64 X 1.32 


256 X 256 X 128 


3.0 X 10- 


-4 


8.0 X 10- 


-5 


1 


No 


0.0 


FDA 


50 



Note. — Col. (1): Name of run. Col. (2): Box size in units of scale heights. Col. (3): Grid resolution. Col. (4): Viscosity coefficient. Col. (5): 
Resistivity coefficient. Col. (6): Dissipation order. Col. (7): Stratification. Col. (8): A-Iean imposed vertical field. Col. (9): Shear advection scheme. 
Col. (10): Total run time in orbits (T^^b = 27rJ7~lj. 



hype rviscosity conserves momentum (jJohansen fc Klahrl 
[2001 . but the energy dissipation is not guaranteed to be 
positive by mathematical construction, nor is the rate-of- 
strain tensor symmetric. We therefore compared all our 
results with the evolution obtai ned by using a strict hy- 
pervi scosity operator f following IHaugen fc Brandenburg 
120041 ) and found no qualitative differences, so we opted 
to use the simplified hyperviscosity because it is much 
simpler and quicker than the full version. 

Shock viscosity {vsh)'- The stratified runs produce 
shocks high up in the atmosphere of the disk. In or- 
der to dissipate sufficient energy in the shocks, without 
rendering the whole simulation domain strongly viscous, 
we have added an additional shock viscosity term. The 
shock viscosity coefficient is obtained in a von Neumann- 
Richtmyer fashion by taking the negative part of the di- 
vergence of u, then taking the maximum over three grid 
cells in each direction, and finally averaging over three 
grid cells, to obtain 

i^sh = Csh(max[-V • u]+){5x)'^ . (8) 

We set the shock viscosity coefficient to Csh = 1-0 in 
the stratified runs to dissipate energy in shocks forming 
above two scale heights from the mid-plane of the disk. 
The shock viscosity Vsh decreases rapidly with increasing 
resolution, so convergence tests also probe the effect of 
decreasing the shock viscosity term. 

2.1.2. Resistivity 

For the resistivity term we consider either normal 
resistivity or hyperresistivity, 

/^^^nV^A + T^aV^A. (9) 

As for the hyperviscosity we use a simplified version for 
V^, but checked that the strict high order Laplacian 
V^V^V^ gives qualitatively and quantitatively similar 
results. 

2.1.3. Diffusion 



Mass diffusion consists of hyperdiffusion and shock dif- 
fusion 

.fo=i?3VV + DshVV + VD,h- Vp, (10) 

where Dgh ~ i^sh is defined in equation ([8|). Mass dif- 
fusion turned out to be necessary in order to damp out 
spurious small-scale modes that are left behind due to 
the dispersion error in the finite differencing of the ad- 
vection term in the continuity equation. We also add a 
bit of extra mass diffusion in the stratified simulations 
to stabilize the shocks forming above two scale heights 
from the mid-plane of the disk (setting the shock diffu- 
sion coefficient Dgh = 1.0 in Eq. [TO]). 

2.2. Boundary Conditions 

We use the usual shearing box boundary condi- 
tions in the radial and azimuthal directions (shear- 
periodic/periodic). At the upper and lower boundaries 
we apply periodic boundary conditions, conserving the 
average fiux of all components of the magnetic field. 
This turned out to be numerically less challenging, and 
we checked that perfect conductor boundary conditions 
(dBx/dz — dBy/dz = B, ^ 0) and vertical field bound- 
ary conditions {Bx = By ^ dB^/dz = 0) for the mag- 
netic field give results that are qualitatively and quanti- 
tatively similar to what we see with periodic boundary 
conditions. 

2.3. Simulation Parameters 

We run simulations for various box sizes, resolu- 
tions, dissipation parameters, and imposed magnetic 
field strengths. The simulation parameters are given 
in Table [TJ The main simulations are the stratified S, 
M, L; and H runs. These all share a vertical extent of 
Lz = 5.28H, but the radial/azimuthal extent of the box 
is varied in steps of two, from = Ly = 1.32 in the 
smallest box^ (S) to = Ly = 10.56 in the largest box 

^ The choice of L = 1.32 as the basic box size unit is based on 
the fact that = 1.32 approximately marks the transition from 
subsonic to supersonic Keplerian shear flow. 
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TABLE 2 

Turbulence Properties 



Run 


























(pUxUy) 




{-BxBy) 




a 


S 


1.9 X 10~ 


-3 


1.7 X 10" 


-3 


1.5 X 10" 


-3 


1.2 X 10" 


_3 


9.1 X 10- 


-3 


5.3 X IQ- 


-4 


1.1 X 10" 


-3 


4.3 X 10" 


-3 


0.0036 


M 


3.9 X 10~ 


-3 


5.1 X 10" 


-3 


2.7 X 10" 


-3 


4.2 X 10" 


-3 


3.6 X 10" 


-2 


2.1 X 10" 


-3 


2.6 X 10" 


-3 


1.2 X 10" 


-2 


U.0U97 


L 


4.0 X 10~ 


-3 


4.5 X 10" 


-3 


2.4 X 10" 


-3 


3.7 X 10" 


-3 


2.6 X 10" 


-2 


2.0 X 10" 


-3 


2.5 X 10" 


-3 


1.1 X 10" 


-2 


0.0089 


H 


4.5 X 10~ 


-3 


4.2 X 10" 


-3 


2.4 X 10" 


-3 


3.5 X 10" 


-3 


2.3 X 10" 


-2 


1.8 X 10" 


-3 


2.6 X 10" 


-3 


1.0 X 10" 


-2 


0.0085 


M_r2 


2.1 X 10~ 


-3 


2.8 X 10" 


-3 


1.4 X 10" 


-3 


2.3 X 10" 


-3 


1.5 X 10" 


-2 


1.2 X 10" 


-3 


1.3 X 10" 


-3 


6.5 X 10" 


-3 


0.0052 


L_nogz 


3.2 X 10" 


3 


2.5 X 10- 


3 


1.4 X 10- 


3 


1.9 X 10- 


-3 


1.3 X 10" 


'2 


7.1 X 10" 


-4 


1.8 X 10" 


'3 


6.9 X 10- 


-3 


0.0058 


M_Bz_r0.5 


5.8 X 10- 


•J 


7.6 X 10- 


3 


4.0 X 10- 


3 


6.1 X 10- 


-3 


6.6 X 10" 


-'1 


3.3 X 10" 


-3 


3.8 X 10- 


-3 


1.8 X 10- 


-2 


0.0148 


M.Bz 


5.9 X 10- 


-3 


1.0 X 10" 


-2 


4.7 X 10- 


-3 


1.1 X 10- 


-2 


7.2 X 10- 


-2 


6.8 X 10" 


-3 


4.1 X 10" 


-3 


2.5 X 10- 


-2 


0.0196 


M_Bz_r2 


5.7 X 10- 


-3 


1.2 X 10" 


-2 


5.0 X 10- 


-3 


1.2 X 10- 


-2 


6.7 X 10- 


-2 


7.8 X 10" 


-3 


4.0 X 10" 


-3 


2.7 X 10- 


-2 


0.0207 


MS 


2.9 X 10- 


3 


2.0 X 10- 


3 


1.2 X 10- 


3 


1.4 X 10- 


-3 


1.0 X 10- 


-'2 


4.6 X 10" 


-4 


1.5 X 10- 


-3 


5.2 X 10- 


-3 


0.0045 


MSj-2 


1.7 X 10- 


3 


1.2 X 10- 


-3 


7.2 X 10- 


-4 


9.1 X 10- 


-4 


5.7 X 10" 


'3 


3.4 X 10" 


-4 


8.4 X 10- 


-4 


3.2 X 10- 


-3 


0.0027 


MSj-2Jix 


3.4 X 10" 


■3 


2.1 X 10- 


-3 


1.7 X 10- 


-3 


1.4 X 10" 


-3 


9.1 X 10" 


-3 


6.2 X 10" 


-4 


1.5 X 10" 


-3 


5.3 X 10" 


-3 


0.0046 


SSj-4_nul 


6.5 X 10- 


-4 


6.1 X 10" 


-4 


4.2 X 10- 


-4 


5.1 X 10" 


-4 


6.3 X 10" 


'3 


1.8 X 10" 


-4 


4.0 X 10" 


-4 


2.8 X 10- 


-3 


0.0021 


SSj-8.nul 


1.0 X 10- 


3 


9.9 X 10" 


-4 


7.1 X 10- 


-4 


8.7 X 10- 


-4 


8.9 X 10- 


3 


3.3 X 10" 


-4 


6.6 X 10" 


-4 


4.2 X 10- 


-3 


0.0033 


MSj-4_nul 


1.9 X 10- 


3 


1.8 X 10" 


-3 


8.7 X 10- 


-4 


1.5 X 10- 


-3 


1.4 X 10- 


-2 


5.1 X 10" 


-4 


1.2 X 10" 


-3 


6.8 X 10- 


-3 


0.0053 



KOTE. — Col. (1): Name of run. Col. (2)-(4): Kinetic energy. Col. (5)-(7): Magnetic energy. Col. (8): Reynolds stress. Col. (9): Maxwell stress. Col. (10) 
o-valuc. Energies and stresses have been normalized with the rrrean thermal pressure in the box, (P) = c^(p). 



(H). For the M run we perform a convergence test by 
doubling the resolution (M_r2), and we also present an 
unstratified version of run L (L_nogz), ignoring the — ^?^z 
term in equation ([T]). The two largest box sizes (runs L. 
L_nogz and H) are evolved using the interpolated Ke- 
plerian shear scheme SAFI, described in Appendix [b1 
This is done in order to relax the relatively tight time- 
step constraint of the Keplerian advection in these large 
boxes and at the same time remove the numerical diffu- 
sivity associated with the Keplerian advection. 

The next group of simulations in Table [1] runs 
M_Bz_r0.5, M_Bz, M_Bz_r2, have a weak vertical mag- 
netic field imposed on the gas, to keep the turbulent ac- 
tivity relatively constant when increasing the resolution. 
This makes it possible to better tell what effect resolu- 
tion has on the zonal flows. In the unstratified runs MS, 
MS_r2 and MS_r2Jix, on the other hand, turbulent ac- 
tivity is maintained when keeping the hyperdissipation 
parameters constant as the resolution is increased. 

The final set of simulations, runs SS_r4_nul, SS_r8_nul, 
MS_r4_nul, apply more natural second derivative viscos- 
ity and resistivity operators to dissipate energy. A rela- 
tively high base resolution is needed here, as otherwise 
there are not enough small scales for the second order 
operators to dissipate kinetic and magnetic energy on. 
Strong zonal flows appear also with second derivative 
dissipation terms, and we see relatively good convergence 
in the zonal flow structure up to our highest resolution 
of around 200 grid points per scale height (eight times 
higher than the S - H runs). 

2.4. Units and Initial Conditions 

Throughout the paper we will state length in units of 
the scale height H, velocity in units of the isothermal 
sound speed Cg, time in units of the inverse Keplerian 
frequency f2~^ or the orbital period Td-b = 27ri7~^, den- 
sity in units of the mid-plane density po and magnetic 
field in units of Cs(AtoPo)^^^j related to the square root of 
the thermal pressure. Energies and stresses are always 
normalized by the mean thermal pressure in the box. 

In the stratified runs we set the initial condition for 



the density through an isothermal hydrostatic equilib- 
rium with scale height H = Cg/f^, whereas the density 
is simply initialized to be constant in the non-stratified 
runs. The toroidal component of the magnetic vec- 
tor potential is set as a moderate scale wave Ay = 
Aocos{kj;x) cos{kyy) cos{kzz) of amphtude Aq = 0.04 
and wave number = ky = kz = 4.76i?~^, yielding 
an initial plasma (3 of around 50. We perturb the gas 
velocity components with random fluctuations of ampli- 
tude Su ~ 10~'^Cs. 

3. TURBULENCE PROPERTIES IN LARGE SHEARING 
BOXES 

The measured properties of the turbulent velocity and 
magnetic fields are written in Table [2] All quantities 
have been normalized by the average pressure in the box, 
(P) = c2(p), for better comparison between stratified 
and non-stratified runs. The last column of Table p]ind i- 
cates the a- value of the flow ijShakura fc Sunvaevlll973f) , 
defined through the Reynolds and Maxwell stress as 

2 {{pu:,uy - B^By)) 

3 (P) ■ ^ ' 

The time evolution of kinetic energy, magnetic energy 
and turbulent stresses for the standard runs S, M and 
L is plotted in Fig. [TJ During the first orbit there is 
a sharp increase in both kinetic and magnetic energies, 
followed by a rapid drop. After the drop the energies 
climb slowly towards their saturation values, which are 
reached after around 20 orbits. In the saturated state 
there is more than a doubling of the energies and stresses 
when increasing the box size from = Ly = 1.32 (S) 
to Lx = Ly = 2.64 (M), but little difference between 
run M and the twice as large run L. The vertical extent 
of the box is fixed at = 5.28 in all the simulations. 
The kinetic energy is distributed relatively isotropically 
in the three components in the small box. But a dis- 
tinct split into a high radial/azimuthal contribution and 
a low vertical contribution is seen when going to larger 
computational domains. This split indicates a transition 
to scales dominated by the Coriolis force, as geostrophic 
modes excited by the turbulence involve only in-plane 



Johanscn, Youdin, & Klahr 



10" 




10"' 



1.32 X 1.32x5.28 



2.64x2.64x5.28 



5.28 x5.28x5.28 



<lpul>l<P> 
<[pi?:>l<P> 




10" 



20 40 60 80 100 20 40 60 80 100 20 40 60 80 100 




20 40 60 80 100 20 40 60 80 100 20 40 60 80 100 



tn... 



10" 



10" 



10"^ 



<PM,M,,>/<P> 

<-BB,>l<P> 




r= 100-200 



10"^ 



20 



40 



60 



80 



100 20 



40 



60 



80 



100 20 



40 



60 



80 



tlT 



100 



Fig. 1. — Kinetic energy (top row), magnetic energy (middle row) and turbulent stresses (bottom row) as a function of time measured 
in orbits, shown for three different box sizes along the columns (standard runs S, M and L). All quantities have been normalized by the 
mean thermal pressure in the box. Energies and stresses increase by more than a factor two when doubling the box dimensions from 



1.32_ff to Lx 



2.64//, but another doubling makes no significant difference for the statistical properties of the turbulence. 



The dotted line in the lower, right panel indicates the further evolution of the Maxwell stress from t 
that the turbulent stresses in run L remain fairly constant after t = 30Torb- 



lOOTorb to t = 200Torb, showing 



gas motion. 

The magnetic field is dominated by the azimuthal com- 
ponent for all domain sizes, as expected since Keplcrian 
shear stretches radial field and creates azimuthal field of 
order \By\ ~ (icohAshoar)|Sx|, where <coh is the coher- 
ence time of the radial field and ishoar = (2/3)f2~^ is the 
shear time-scale. The ratio rms(i3y)/rms(i?x) ~ 2.5-3.0 
for runs S-H gives a (scale averaged) coherence time-scale 
of ^coh ~ 2J7~^. Based on the long correlation time of 
large scale density and velocity structures , one might 
expect an even greater dominance of By, especially in the 
larger boxes. However, the large scale zonal flow is built 
up from small, uncorrelated kicks by the Lorentz force, 
and a relatively short correlation time of the magnetic 
field fits well into this random walk picture. 
Maxwell and Reynolds stresses increase from 



smallest box to (pu^Uy) w 0.0025 and {—BxBy) w 0.01 
in the largest box, yielding an increase in a-value 
from a w 0.004 to a « 0.009. This increase is of 
interest because global simulations of the MRI tend to 
yield hi gher g-values than tr aditional narrow shearin; 



boxes (IHawlev et al 
lArlt fc Riidiged 12001 



0.001 and (-B^B., 



0.004 in the 



19951 iBrandenburg et all [l99 
Fromang fc NelsonI \200¥ . Our 
simulations suggest a transition towards the global 
value with increasing size of the shearing box. An 
additional explanation for the measured high stresses in 
global simulations may be the relatively low resolution 
per scale height possible in these simulations. The 
saturated state of non-stratified zero net vertical flux 
magnetorotational turbulence ha s been shown to decline 
with increasing resolution (jFromang fc Papaloizoul 
|2007[) . due to the change in numerical dissipation 
as the resolution increases. The turbulent transport 
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Fig. 2. — The kinetic (blue/dark) and magnetic (yellow/gray) 
energy spectrum, summed over shells of constant wavenumber k = 
\k\. The kinetic energy dominates at the largest scales in the small 
box (top panel), whereas there is equipartition between kinetic and 
magnetic energy in a big range of smaller scales. In the two larger 
boxes (middle and bottom panels) the magnetic energy comes to 
dominate scales larger than approximately two scale heights. At 
the same time the equipartition at small scales is broken in favor 
of a dominance of kinetic energy. 

coefficients in run M_r2, with twice as high resolution 
as run M, are approximately 50% lower than in run 
M (see Table [2]). Such a linear scaling of transport 
coefficients with resolu t ion ex tends the effect found by 
iFromang fc Papaloizoul ()2007f ) to stratified boxes. 

The kinetic and magnetic energy spectrum arising in 
boxes of different sizes is shown in Fig. [21 summed over 
wavenumber intervals I A; I g [fccen— fccGn+fco/2], with 
fccon taking integer values times the smallest wavenumber 
fco. Wc have only included isotropic scales in the summa- 
tion. Thus we exclude (a) long vertical wavelengths that 
exceed the radial or azimuthal extent of the smaller boxes 
in runs S and M, and (b) the anisotropic high wavenum- 
bers which only exist in the "corners" of fc-space. We 
averaged the energy spectra over 71 snapshots taken be- 
tween 30 and 100 orbits. 

For the two large runs (M and L) the kinetic en- 
ergy (blue/dark) line dominates over magnetic energy 



(ycllow/gray) at scales below approximately 1 ... 2 scale 
heights. At larger scales (smaller k) the kinetic energy 
falls off gradually, whereas the magnetic energy keeps ris- 
ing until finally peaking near the largest scales of the box. 
This transition from kinetically to magnetically domi- 
nated flow can only be seen when considering a suffi- 
ciently large box size. Interestingly the magnetic energy 
is comparable to kinetic energy near the dissipative sub- 
range in the smallest box, whereas the two energies move 
apart when the box size is increased. 

In Fig. [3] we show for run L the spectral power of all 
dynamical components [u^, Uy, Uz, p, Bx, By, B^) as 
a function of radial wavenumber (cc-axis), the smallest 
two azimuthal wavenumbers (columns) and the smallest 
three vertical wavenumbers (colored lines) in the first two 
columns. The third column, on the other hand, shows 
the spectral power as a function of vertical wavenumber, 
but for different, fixed values of the radial wavenumber. 
Looking first at the magnetic field (three bottom rows of 
Fig. [3]) one sees that the radial and azimuthal magnetic 
energy peaks at the largest axisymmetric scales of the 
box [ky = 0), whereas the non-axisymmetric modes have 
far less power. However, the large scale purely vertical 
modes shown in the third column contain approximately 
a factor four times more magnetic energy than the radial 
modes. The vertical component of the magnetic field has 
a rather flat dependence on all directions and show little 
large scale excitation. 

The azimuthal velocity field component and the gas 
density also exhibit very high power at large axisymmet- 
ric scales; we return to those zonal flows in 21 

3.1. Magnetic Fields on Large Vertical Scales 

The large scale power of the z-dependent radial and 
azimuthal magnetic field components, seen in the last 
column of Fig.[3l is illustrated in Fig.[4l In this figure we 
show the azimuthal magnetic field, averaged over the x 
and y directions, as a function of time t and height over 
the mid-plane z. Magnetic field structures propagate 
away from the mid-plane such that the field at a given 
height over the mid-plane flips signs at regular inter- 
va ls. Similar buoyancy wave s are seen in the simula tions 
bvlBrandenburg et all (flOOl . iMiUer fc Stond (|2000t) and 
iHirose et al.l (|2006f ). The period of the sign flip is around 
13 orbits, in agreement with the turbulent diffusion time- 
scale = H^/rjt, where T]t ^ 0.01 is the turbulent re- 
sistivity. This time is much longer than the correlation 
time of less than one orbit deduced in the previous sec- 
tion from the creation of azimuthal field from the radial 
field by the Keplerian stretching term. 

The vertically propagating radial and vertical field may 
be the effect of a dynamo driv e n by rotation , strati- 
fication and shear (jParked ITOSSl : iMoffattI [iJTSh . How- 
ever, by correlating the mean field with the fluctuat- 
ing electromotive f o rce at both sides of the mid-plane, 
[Brandenburg et all ()1995D concluded that the dynamo 
coefficient has the opposite sign of what is expected from 
mean field dynamo theory. 

The non-stratified simulation L_nogz displays verti- 
cally dependent magnetic field structures propagating 
from and to the mid-plane. Without stratification there 
is no systematic helicity separation over the mid-plane, 
although patches of positive and negative helicity ap- 
pear with a correlation time of a few orbits. The large- 
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(two first columns) and vertical wavenumber (last column). Colors indicate different wavenumbers in either the vertical 
direction (two first columns), or in the radial direction (last column), while the azimuthal wavenumber ky is indicated at the top of each 
column. There is a clear tendency for Ux, Uy, p, Bx and By to peak at the largest radial scale of the box. However the last column shows 
that most large scale magnetic energy resides in modes with vertical variation. The variation of p with z is mainly due to the stratification. 



scale magnetic field activity m ay come around as a result 
of an incoherent a-oj-dynamo ((Vishniac fc Brandenburg 
I1997I) . tliriving on tlie fluctuating helicity, or be con- 
nected tojnBgnfitoroteti£nal instability in the azimuthal 
fields (E esur fc Ogilvid[200l . 

4. ZONAL FLOWS 

In Fig. [5] we show the gas density at the sides of the 
simulation box at t = 72 Torb for run L (box size of 
5.28i? X 5.28i7 x 5.28iJ). The turbulent structure of 
the gas density is clearly seen superimposed on the main 



stratification in the left panel. The right panel shows the 
column density along each direction Ui = J p(x, y, z)Axi^ 
divided by the directional column density in the initial 
stratification. There is a dominant column density struc- 
ture at the largest radial scale of the box, with very little 
dependence in the azimuthal and vertical directions'^. 
The space-time plots of the gas density of runs S, M 

^ The marked drop in column density near the upper and lower 
boundary planes is likely due to the periodic boundary conditions 
that precludes hydrostatic equilibrium near the boundaries by forc- 
ing the density derivative to be zero. 
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Fig. 4. — The azimuthal magnetic field By of run L as a function 
of time and height over the mid-plane. Field structures continu- 
ously propagate away from the mid-plane, causing the azimuthal 
field at a given height over the mid-plane to flip signs with a period 
of approximately 13 orbits. 

and L shown in Fig. [6] confirm the dominance of the 
kx/ko ~ 1 mode. Here we have averaged over the vertical 
and azimuthal directions and show the gas density as a 
function of radial coordinate x and time t measured in 
orbits. Both the density amplitude and the correlation 
time of the large scale structures increases significantly 
when making the box size larger. Run H is not shown in 
Fig. [5]- the long correlation time of the density structures 
in this huge box would require several hundred orbits in- 
tegration to get good statistics. 

Axisymmetric density modes have their power peak at 
the largest radial scale of the box. The azimuthal veloc- 
ity has a similar peak at the largest axisymmetric modes 
(see Fig. Fig. [7] compares the large-scale behavior 
of the density fluctuations to other physical quantities. 
Quantities are radially shifted so that the k^ = ko den- 
sity bump peaks at the plot center, and then are time 
averaged. This filtering causes the plotted quantities (es- 
pecially the gas density) to appear sinusoidal on the box 
scale, since higher order modes do not have a consistent 
phase relationship with the largest mode. 

Fig. [7] shows an almost exact — 7r/2 phase difference 
between the radial peaks of the density and azimuthal 
velocity. This follows from near geostrophic balance be- 
tween Coriolis forces and pressure, 

2pQf2uyKdPldx, (12) 

which gives super-Keplerian (sub-Keplerian) motions ra- 
dially interior (exterior) to the density peak. In detail 
though the radial force balance includes Lorentz forces. 
The magnetic pressure and Maxwell stress perturbations 
in Fig. U] are almost perfectly anti-correlated with the 
thermal pressure (the trend for thermal and magnetic 
pressures to misalign comes naturally out of the simpli- 
fied zonal flow model presented in fj5]). This supports 
recent findings by T. Sano (personal communication. 



2008). 

Azimuthally and vertically varying modes have been 
averaged out of the density structures that appear in 
Fig. El But Fig. [3] shows that there is a clear trend 
for the density power to fall quickly with azimuthal 
wavcnumber, likely because the Keplerian shear drives 
non-axisymmetric spectral power towards larger and 
larger k^- The decreased spectral power of density at 
vertically varying modes can be seen as an effect of the 
Taylor-Proudman theorem. There is no geostrophic bal- 
ance available for modes with vertical variation other 
than the main stratification, so pressure and turbulent 
diffusion will even out any vertical variation faster than 
the magnctorotational instability injects energy to those 
modes. 

In Table [3] we show the measured values for the Fourier 
amplitudes of azimuthal velocity, density and Maxwell 
stress at the largest radial scale of the box, the power law 
index for the two largest radial scales, and the correlation 
time of the density structures. The power law index is 
found from 

^_ ln[|/|(2fco,0,0)]-ln[|/|(fco,0,0)] 

^~ hr2 ■ ^^-^^ 

We have normalized the density amplitude by the mean 
density in the box to get effectively the column den- 
sity perturbation relative to the mean column density, 
/5(A:o,0,0)/(/9) = S{ko,0,Q)/{S). The correlation time 
has been calculated by taking the value of p, averaged 
over y and 2, at a given time t and measuring for each 
spatial point the time it takes for the density at that 
point to change by a value corresponding to the stan- 
dard deviation of the gas density. The result has been 
averaged over many closely spaced times (for which the 
turbulence is saturated and the correlation time does not 
extend to the final time of the simulation) and multiplied 
by two to represent the full temporal extent of the cor- 
related structures. The resulting correlation times are in 
good agreement with the life-time of the structures seen 
Fig. [6l and the measurements confirm the very long cor- 
relation times seen in simulation L. For the even larger 
box, simulation H, the correlation time is approximately 
50 orbits, comparable to the total simulation time. 

Interestingly the amplitude of the large scale azimuthal 
velocity is relatively independent of box size, while the 
density amplitude grows approximately proportional to 
the radial extent of the box"^. There is an almost per- 
fect geostrophic connection between the amplitude of the 
zonal flow and the normalized density amplitude, follow- 
ing the expected \uy\ ~ {l/2)ko\p\. An exception to this 
is the mean vertical field runs (M_Bz_*) where the az- 
imuthal velocity is higher than needed for geostrophic 
balance. We give a possible explanation for this discrep- 
ancy in i j5.41 based on the relatively short correlation 
time of the turbulence in the imposed field runs. 

The power law index of azimuthal velocity and density 
also follows the geostrophic expectation of Cp = C^^ 
although both and Cp decrease (towards zero) with 
increasing box size. This decrease could indicate some 

* Parceval's theorem would imply that the power at a given 
Fourier mode is inversely proportional to the number of grid points, 
but for the case of zonal flows, the power peaks so strongly at the 
largest radial mode that this reduction does not apply. 
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Fig. 5. — The gas density at the sides of the simulation box (left panel) and the directional column density (right panel, normalized by 
the column density of the original stratification) of run L at time t = 72rorb ■ The dominant column density structure is vertically constant 
and axisymmetric and has variation mainly at the largest radial scale of the box. 
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Fig. 7. — Phase-shifted and time-averaged perturbations to the 
gas density, azimuthal velocity, Maxwell stress (multiplied by 10 for 
visibility) and magnetic pressure (also multiplied by 10) for run L, 
as a function of the radial phase <!>. The —it/ 2 phase shift of Uy 
is a consequence of geostrophic balance on large scales. There is 
a clear anti-correlation of gas pressure (proportional to density) 
with Maxwell stress and magnetic pressure. Perturbations are 
azimuthally and vertically averaged for each time-step. Prior to 
time-averaging, the quantities are shifted in radius so that the low- 
est order radial perturbation (with kx = ko = 2tc/Lx) to the gas 
density is centered in the plot. 

motion towards convergence of the zonal flow and pres- 
sure bumps. Global simulations will be needed to eventu- 
ally find such convergence, since curvature effects could 
play an important role for simulation domains of more 
than ten scale heights. 

The presence of large scale radial variation in Maxwell 
stresses is clear from Table [3l The amplitude is typically 
around 10-15% of the mean Maxwell stress, with a sharp 
decrease towards shorter wavelengths {(s^By ~ — !)■ 
The imposed field simulations (M_Bz_*) even indicate 
an increasingly top heavy Maxwell stress with increas- 
ing resolution. However, the turbulence activity in that 
set of runs also increases by 50% as the resolution is in- 
creased (see Tabled, while the amplitude of the large 
scale Maxwell stress increases somewhat slower. 

4.1. Large-Scale Balance 

To get an understanding of the launching mechanism 
for zonal flows, we plot in Fig. [8] the energy balance of 
all vertically constant and axisymmetric modes. We have 
taken each term Qi in the dynamical equation of variable 
f,f~ "^iQi, and calculated its contribution to the power 
of all the purely radial modes as 



d\fl 
dt 



= ^2Re(.r4:). 



(14) 



Here hats denotes Fourier transforms and stars com- 
plex conjugation. Wc add the power contribution from 
both the negative and the positive wavenumber in Fig.O 
The second panel of Fig. [5] indicates that the zonal flow 
[uy{kx)\ is excited primarily by the Lorentz force (indi- 
cated by the positive value of the "Lor" curve). A large 



scale radial variation in the Maxwell stress (see Table [3] 
and Q exerts a significant azimuthal tension on the gas. 
As the zonal flow is launched by the magnetic tension, 
the Coriolis force creates a slight radial component of the 
velocity fleld. The ensuing convergence of the gas flow 
creates the radial pressure bumps seen in Fig.[B]in perfect 
geostrophic balance between Coriolis force and pressure 
gradient force. 

Fig. [9] shows that there is some correlation between 
the phase of the zonal flow and the phase of magnetic 
tension associated with the Maxwell stress, although the 
latter varies on a much shorter time-scale than the zonal 
flow. The magnetic tension is only an order of magnitude 
lower in amplitude than the saturated zonal flow, so a 
correlation time of the magnetic field of a few orbits is 
enough to explain the build up of large scale flow. 

5. SIMPLIFIED MODEL OF ZONAL FLOWS 

The zonal flows in our simulations have a very long 
correlation time, on the order of many tens of orbits. We 
can understand both the amplitude and the life-time of 
the zonal flows by adopting a simplified dynamical model 
which is forced by the azimuthal magnetic tension. Con- 
sider the following set of equations for the evolution of 
the radial velocity Ux, azimuthal velocity Uy and density 
p, all for the large scale k = (fco, 0, 0) mode, 



Q = 2nuy -ikf)p, 



Pa 



duy 1 . 

— - = ilux - 

dt 2 

'^P -i, - 

dt 



(15) 
(16) 
(17) 



Equation ([T5|) expresses radial geostrophic balance^. We 
ignore for simplicity the role of Lorentz forces in this bal- 
ance, evident from Fig. [5] to be an order unity correction. 
The azimuthal force balance in equation (|16p describes 
forcing by the azimuthal tension T, which is partly de- 
fiected into radial motion, but also directly accelerates 
the zonal flow in Uy. Mass continuity in equation (jl7p 
includes a damping term with the timescale Tmix. This 
represents non-linear advection terms, seen to be signifi- 
cant in Fig. [5] (bottom left panel). The same figure (top 
left and top center panels) shows that non-linear advec- 
tion makes a negligible contribution to the momentum 
equations. 

We can combine the above equations into a single evo- 
lution equation for the density, 



dp 
dt 



1 



1 



F 



(18) 



where F = ^2ik^pQT / Q is the forcing term. The prefac- 
tor Cfc = (1 -I- kgH^)^^ is a pressure correction for small- 
scale modes that both decreases the amplitude of the 
forcing and increases the effective damping time. The 
density bump grows in exact antiphase with the stress 
bump S, equivale nt to the = const relation of steady 
global accretion ()Pringldll98lD . This follows from writ- 

^ Dropping dux/dt filters high frequency density waves which 
are not relevant for the longer-lived structures considered here. 
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TABLE 3 
Zonal Flow Properties 



Run 






\fi,,\(kn 01 


\h\(kn 01 




ffTB^Kkn 01 




sp 


By 


f'corr 


S 


4.3 X 10" 


3 


1.9 X 10"2 


8.1 X 10- 


3 


4.5 X 10—* 


-1.6 


-2.4 


-0.9 


8.3 


M 


1.2 X 10- 


2 


2.3 X 10-2 


1.9 X 10- 


2 


1.2 X 10-3 


-1.0 


-1.8 


-0.8 


8.5 


L 


1.1 X 10" 


2 


1.7 X 10-2 


3.3 X 10- 


2 


1.3 X 10-3 


-1.0 


-2.3 


-1.4 


19.6 


H 


1.0 X 10" 


2 


1.5 X 10-2 


5.2 X 10- 


2 


1.1 X 10-3 


-0.2 


-1.1 


-0.9 


50.5 


M_r2 


6.5 X 10- 


3 


1.6 X 10-2 


1.3 X 10- 


2 


5.5 X 10-'' 


-0.9 


-2.2 


-1.0 


12.6 


Ljiogz 


6.9 X 10" 


3 


1.2 X 10-2 


2.1 X 10- 


2 


5.6 X 10-'' 


-0.8 


-1.7 


-1.0 


18.8 


M_Bz_r0.5 


1.8 X 10" 


2 


4.5 X 10"^ 


2.4 X 10- 


2 


2.1 X 10-3 


-1.1 


-1.5 


-0.6 


4.6 


M.Bz 


2.5 X 10" 


2 


3.9 X 10-2 


2.4 X 10- 


2 


2.4 X 10-3 


-0.6 


-1.9 


-0.9 


5.4 


M_Bz_r2 


2.7 X 10" 


2 


3.7 X 10-2 


1.8 X 10- 


2 


2.5 X 10-3 


-0.5 


-1.5 


-1.2 


3.2 


MS 


5.2 X 10" 


3 


2.7 X 10-^ 


2.3 X 10- 


2 


9.6 X 10-" 


-1.1 


-1.5 


-0.9 


9.9 


MSj-2 


3.2 X 10" 


3 


1.5 X 10-2 


1.3 X 10- 


2 


4.0 X 10-" 


-0.8 


-1.3 


-0.8 


11.3 


MSj-2Jix 


5.3 X 10" 


3 


2.9 X 10-2 


2.5 X 10- 


2 


7.7 X 10-" 


-0.9 


-1.4 


-0.5 


15.7 


SS_r4_nul 


2.8 X 10" 


3 


2.0 X 10-^ 


8.5 X 10- 


3 


3.9 X 10-* 


-1.3 


-2.1 


-0.8 


6.5 


SS.r8.nul 


4.2 X 10" 


3 


2.2 X 10-2 


9.4 X 10- 


3 


5.4 X 10-" 


-1.5 


-2.3 


-0.9 


7.6 


MSj-4_nul 


6.8 X 10" 


3 


3.3 X 10-2 


2.8 X 10- 


2 


9.3 X 10-" 


-1.7 


-2.4 


-1.0 


9.4 



Note. — Col. (1): Name of run. Col. (2): Maxwell stress, normalized by mean pressure. 
Col. (3)-(5): Fourier amplitude of azimuthal velocity, density normalized by mean density in the 
box, and pressure-normalized Maxwell stress, respectively, at largest radial mode. Col. (6)-(8): 
Power law index of azimuthal velocity, density, and Maxwell stress, respectively, calculated from 
two largest radial modes. Col. (9): Correlation time, in orbits T = 2nn~^, of largest radial 
density mode. 



d/dt[k^ul(k^,0,0)] 
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-5- 10-' 



1 



10 



kJkr, 



d/dt[k^(k^,0)] 



1 1 1 — I — ■ ■ ■ ■ 1 ■ ■ ■ — ■ — ■ 

Adv — 

V Cor — 
■\, ^ Lor ■ 
\ \ ^ Pre - - 


1 1 1 1 ■— I— i-i-i 1 1 1 1 r 

Adv — 
Cor — : 
Lor — 


/ 





10 

kJkn 




10 



kJkr. 



hlQ- 



-hlO' 



d/dt[k^p\k^,Om d/dt[k^Bl(k^,Om 





Adv — 
Com — 

-vvi^/t 








Adv — . 
Com — 
Str 



10 



kJk 







10 

kjkr\ 



dldt[k^l{k^,Qm 



Adv — 
Com — 
■ Str — 









10 



kJkr. 



1-6 



MO 
5- 10"' 



-5- 10"' 
-MO"^ 



Fig. 8. — Power budget for all axisymmctric and vertically constant modes in run L. The curves show the power contribution from 
individual terms in the dynamical equations ["Adv"=advection, "Cor"=Coriolis force, "Lor"=Lorentz force, "Pre"=pressurc gradient, 
"Com"=compression, "Str" =stretching] . The data has been averaged over equally spaced snapshots from t = SOT^rb to t=100Torb and 
multiplied by kx for better visibility of the small scale power. The energy in the azimuthal velocity field (i.e. the zonal flow) is pumped 
by the Lorentz force (magnetic tension). The resulting convergence in the radial velocity leads to an increase in density amplitude and a 
decrease in azimuthal magnetic field. The compressive increase in gas density is in turn balanced by turbulent diffusion. The non-linear 
advection term in the induction equation for the azimuthal field, on the other hand, adds magnetic energy at the largest scale. 
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Fig. 9. — Phase and amplitude of the large-scale mode 
{kx,ky,kz)/ko = (1,0,0), for the azimuthal component of the ve- 
locity field and the magnetic tension associated with the Maxwell 
stress. There is some correlation between the phases (top panel), 
corroborating that zonal flows are indeed launched by the x- 
dependent Maxwell stress, although the Maxwell stress has vari- 
ation on a much shorter time-scale than the zonal flows. The 
amplitude of the magnetic tension is approximately an order of 
magnitude smaller than the azimuthal velocity (bottom panel). 

ing F = —2kQf2^^S, where the complex stress amplitude 

is defined as = —fiQ^BxBy. 
The straightforward equilibrium solution to equation 

IHI) is pcq ~ Tmix-F- This is however not the correct 
solution to the problem since Fig. [5] clearly demonstrates 
that the magnetic tension forcing T does not maintain a 
consistent phase relationship with Uy (and thus with p) 
during the life-time of a zonal flow. 

Instead equation (fT5|) should be modeled as a stochas- 
tic differential equation that gives the growth of the den- 
sity and azimuthal velocity as a damped random walk. 
We now give a description of the approximate solution^ . 
We assume that the forcing term CkF has a correlation 
time Tfor- The density amplitude changes during each co- 
herent step by |/5stcp| ~ Ck\F\T[oi-, and grows as a random 

^ More formal solution method s exist as a man ifestation of the 
fluctuation dissipation theorem. lYoudin fc Lithw ick (2003) con- 
sider a mathematically similar problem: the stirring of particle ve- 
locities by aerodynamic coupling to turbulent gas. They show that 
the approximate arguments given here are consistent with more 
formal solutions to the stochastic Langevin equation via Fokker- 
Planck analysis. 



walk, p \pstep\y/Nstop ~ IPstoplW'^for- Thc damping 
term sets the lifetime of coherent structures as 

Tcorr ~ Tn,i,/cfc ~ (1 + klH^) / [klO,) , (19) 

where the final step describes the mixing process as dif- 
fusive with a strength Dt. The random walk saturates 
after A^stcp ^ Tcon/Tfor steps, resulting in the equilibrium 
value 

^ = 2\/CfcTforTniix-H"fco— . (20) 

We appeal to the simulations for appropriate values of 
Tmix and Tfor to check the consistency of the model. 
The forcing by magnetic tension can be estimated for 
run L from Fig. [9] to be approximately T ~ 0.001i7cs. 
From the same figure we read the correlation time of 
the forcing to be approximately Tfor ~ 101?^^ (around 
two orbits). We assess the turbulent mixing time-scale 
as a diffusive process with r^ix = l/(fco-Dt) ^ 70n~^. 
Here we use Dt ~ aH'^Q ^ O.OliJ^i? for the turbulent 
diffusion coefficient. 

With the above quantities, equation ((20)) gives an equi- 
librium density amplitude p ^ 0.041 for run L, fairly 
close to the measured value p ^ 0.033. This agreement 
is remarkable for such a simple model and actually a bit 
fortuitous, since Fig. [S] (top left panel) shows that radial 
Lorentz forces (neglected in this toy model) balance a 
significant portion of the pressure gradient, requiring a 
yet higher density perturbation. 

5.1. Life-times of Pressure Bumps and Zonal Flows 

Zonal flow lifetimes are a useful constraint on our un- 
derstanding of large scale dynamics and mixing in MRI 
turbulence. The correlation times in Table [3] can be 
compared to the simple predictions of the random walk 
model in equation pop ^. In large shearing boxes, with 
k^H 1, the lifetime increases with the diffusive mix- 
ing timescale, i.e. icorr oc L^. For smaller boxes, with 
koH ^ 1, the predicted tcorr is independent of box size 
for a fixed diffusion coefficient. 

Indeed, tcorr ~ 8 Torb is nearly constant in the two 
smallest boxes, S and M, in rough agreement with the 
simple model. For runs L and H our estimate gives 
tcorr oc L].'^^ in good agreement with measured value 
of tcorr fx L]:^'^ for these runs. 

5.2. Consistency Checks on Forcing Model 

We now check the consistency of our model further, 
first by computing the expected forcing timescale (the 
least well-measured quantity) in terms of better mea- 
sured quantities. Inverting equation (|20|) gives 

" (^) i ' ^''^ 

where we use i? = Cg = po = 1 and Tcorr = 27ri7tcorr = 
Tmix/cfe. Since these are all known quantities (using 

T oc koBxBy normalized to T = 10~^ in run L), we 
compute tfor = Tfor/(27r) « {0.64,0.62,1.5,9.8} orbits 
for runs {S,M,L,H}. That tfor ^ tcorr supports the basic 
assumption of a random walk model. 

^ Where tcorr = Tcorr /(Stt) since the former is in units of T^ih- 
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As another consistency check, we can measure the tur- 
bulent mixing coefficient, Dt = (Tmixfco)~^ using the re- 
sult from the random walk model that 27rtcorr ~ T^ix/ck- 
We get 

1 + 

Dt = ^^-^ w {2.0,2.2,1.4, 1.2} X 10"^ (22) 

for runs {S,M,L,H}. The L and H runs give reasonable 
values, but the S and M runs have anomalously high mix- 
ing coefficients relative to other indicators of the strength 
of turbulence. However, the random walk model overall 
appears physically consistent and is useful in interpreting 
the simulations. 

5.3. Extrapolation to scales 3> H 

It is interesting to speculate how zonal flows would 
behave in even larger boxes than considered here. This 
is relevant in discs with aspect ratio H/r ^ 0.1, such 
as the di scs around supermassive black holes in active 
galaxies (jFrank et a l. 2002; Goodman 2003). 

Part of the observed growth in density amplitude with 
increasing box size is the corresponding increase of the cu 
factor [eq.[20]. Since Ck tends towards unity with increas- 
ing wavelength, it will not contribute to the increase in 
p in boxes much larger than ten scale heights. Assuming 
further that Tfor and Tmix scale approximately as l/fcg, 
and that T oc fco , we predict that the density amplitude 
will not continue to grow, but will remain relatively con- 
stant on scales A ^ IQH . 

However, this prediction is uncertain, given that we 
do not have a full understanding of how the correlation 
time of the forcing depends on scale. Our analysis in the 
previous section indeed showed that Tfor grew faster than 
1/fco from box L to box H. Also it is not known whether 

BxBy will remain constant at very large scales. Instead 
the inverse cascade of magnetic energy may eventually 
terminate. 

In larger boxes a wider range of scales can contribute to 
the density bumps, because of the long mixing time-scale, 
not just the largest scale that we considered above. This 
"oligarchic" regime may be already have been reached in 
our run H {Lx — 10.56-ff ) where the zonal flow amplitude 
is slightly smaller than in run L, while the zonal flow 
peaks much less strongly at the largest scale {C,uy = — 1.0 
for run L, but is only C,uy = —0.2 for run H). 

5.4. Towards a More Complete Theory 

A more complete theory of large scale dynamics would 
include the self-consistent evolution of magnetic fields. 
For instance, the magnetic pressure grows from the 
Maxwell stress, via the term 

d[{l/2)Bl] 3 , , 

y ^...--QBxBy, (23) 

which readily follows from the induction equation ([5]). 
Thus the magnetic pressure grows in regions of stronger 
(more negative) Maxwell stress. The density, on the 
other hand decreases in regions of high stress. Fig. [7] 
clearly shows the anticorrelation of density (and thus 
thermal pressure) with both Maxwell stresses and mag- 
netic pressure. Though not shown in the plot, actually 



all components of the magnetic pressure are found to be 
lower in the density peak. 

The departure from geostrophic balance for the im- 
posed field runs M_Bz_* likely arises from the relatively 
short correlation time of the more vigorous turbulence in 
these simulations. This can break the validity of equation 
p5|) by launching density waves. From equations p6p 
and p?|) we see that zonal flows can be launched imme- 
diately by the magnetic tension {duy/dt = T), without 
leaking into the Coriolis force {~f2ux/2) that gives ra- 
dial compression. This way a pressure bump does not 
have enough time to react to the zonal flow before the 
magnetic tension shifts position. 

The recent work by IVishniaj ()2008[ ) offers a model 
for the growth of large scale structures in MRI turbu- 
lence, including an induction equation for fleld evolu- 
tion. Other features which differ from our model include 
incompressibility and the neglect of azimuthal accelera- 
tions (while we ignored radial ones). Most important, 
iVishniaj (|2008f ) considers only disturbances with a si- 
nusoidal vertical dependence, which leads to the growth 
of vertically varying rolls. The relation between these 
structures and the height-independent zonal flows that 
emerge from our model deserves further study. 

6. INVERSE CASCADE 

So far we have shown that the large scale variation in 
the Maxwell stress launches slightly compressive zonal 
flows, which in turn enter geostrophic balance with ax- 
isymmetric column density bumps. 

The magnetic energy in the azimuthal field is pumped 
to the largest radial scales of the box by the advection 
term (in the induction equation, see Fig. [8]) . This inverse 
cascade is also evident from Fig. [101 Here we plot the 
power contribution from different terms in the induction 
equation to the magnetic field, summing over shells of 
constant wavenumber and averaging over equally spaced 
snapshots during the saturated state of run L. The ra- 
dial and azimuthal field both get a net positive contri- 
bution from the advection term at large scales, whereas 
the magnetic stretching term (including stretching by the 
background shear) adds energy with a peak at middle 
scales around k/ko ^ 20. The source Bx for Keplerian 
stretching is itself generated by non-Keplerian stretch- 
ing at both moderate and small scales. For the verti- 
cal field component Fig. \W\ shows no sign of large scale 
magnetic energy activity. This is interesting in connec- 
tion with a possible accretion disk dynamo, because the 
lack of large scale activity of B^ may indicate that verti- 
cal field, crucial for the vertical field MRI to operate, is 
only replenished at scales near the dissipative subrange 
of the turbulence, but not on the large scales. Azimuthal 
fields are unstable to a non-axisymmetric magn etorota- 
tional swing instability ()Balbus k HawlevI [1992 ) . which 
has a significant growth rate at short vertical wavelength s 
(jFoglizzo fc Taggedll995t iTerguem fc Papaloizoull 19961) . 
Fig. [To] indicates that the azimuthal field is actively pro- 
cessing energy, through the various terms in the induc- 
tion equation, even at the largest scales of the box. 

Fig. \TU\ also shows that compression is a sink of mag- 
netic energy at the largest scales, for both radial and 
azimuthal field. The zonal flow launched by the radially 
dependent Maxwell stress indeed creates a diverging ve- 
locity field in regions where the magnetic field is high. 
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Fig. 10. — Contribution from different terms of the induction 
equation to tfie magnetic energy at scale k, summed over shells 
of constant k = \k\, including the contribution from the hyperre- 
sistivity ["Res"]. Both the x- and y-components of the magnetic 
field get a positive energy injection at the largest scale from the 
advection term ["Adv"]. Another important energy source is the 
Keplerian stretching term ["Str (K)"]. The regular stretching term 
(["Str"]) injects energy at all scales of Bx and Bz, but is an energy 
sink for By . 



lowering the magnetic energy at the largest scales. This 
way the zonal flows are self-limiting, as they reduce the 
large scale variation in the magnetic field. 

When the magnetic energy in the radial field compo- 
nent cascades up to the largest scales of the box, then 
the Maxwell stress at those scales can grow through the 
negative definite Keplerian stretching term, 



d{B,By)/dt = ...- {3/2)f2Bl 



(24) 



The deepest source of the zonal flows thus lies primarily 
in the inverse cascade of the radial magnetic field energy. 
In Fig. [TT] we show the radial magnetic energy and the 
IVIaxwell stress in run L, averaged over the azimuthal and 
vertical directions and normalized by the instantaneous 
mean. There is an astonishing correlation between the 
two quantities. 





xlH 

B^Bp,t)/<B^B^it)> 




Fig. 11. — The radial magnetic energy (top panel) and the 
Maxwell stress (bottom panel) in run L, both averaged over the 
azimuthal and vertical directions and normalized by the instanta- 
neous mean. There is an almost perfect correlation between the 
two quantities, arising from Keplerian stretching of the radial mag- 
netic field. Energy and stress fluctuations peak at 30% of the mean 
value, but the average fluctuation is only around 10%. The coher- 
ence time is generally a few orbits, with periods of coherence times 
up to ten orbits. The correlation time of the density structures 
shown in Fig. [6] is much longer. 



6.1. Hydrodynamical Inverse Cascade 



9 



The slightly positive energy contribution of the advec- 
tion term to the velocity field (Fig. [5]) represents evidence 
for an additional inverse cascade for the azimuthal veloc- 
ity. Suc h an inverse cascade occurs com monly in 2-D tur- 
bulence (jKraichnanl '1967: 'Boruc '19931 an d in 3-D rotat- 
ing t urbulence (e.g. .Bartello ot a l. 1994; iCambon et al.l 
Il997| ). In the latter case vertically symmetric cyclones 
and anticyclones emerge at scales large enough for Cori- 
olis forces to dominate (i.e. Rossby number Ro ^ 1). 
Cascade of small scale turbulent energy to large scale 
Rossby waves and zonal jets in protoplanetary discs was 
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proposed bv lSheehan et all (|1999f ). However, in our sim- 
ulations the hydrodynamical inverse cascade is found to 
give a more than 10 times smaller contribution to the 
large scale kinetic energy than the variation in Maxwell 
stress and magnetic energy do. 

7. NUMERICAL ISSUES AND VALIDATION TESTS 

Supersonic shear flow in large shearing boxes intro- 
duces a space-dependent numerical dissipation which can 
lead to artifi cial suppression of turbulence near the edges 
of the box (jJohnson et al.l i2008, ) . In this section we 
present various tests to validate the results of the non- 
linear simulations. 

7.1. Results for Various Shearing Box Algorithms 

In Appendix [C] we show variations of runs S-H where 
we use several different ways to solve the shearing box 
equations. The various schemes - FDA (Finite Difference 
Advection), RRC (Random Revolution Center). SRD 
(Systematic Radial Displacement), and SAFI (Shear Ad- 
vection by Fourier Interpolation) ~ are discussed in detail 
in Appendix [C] The statistical properties of the zonal 
flows and the pressure bumps are given in Table HI The 
large scale variation in the Maxwell stress remains at 
around 10% of the mean Maxwell stress, independent of 
the shearing box algorithm. The amplitudes of azimuthal 
velocity and density, on the other hand, are clearly larger 
when applying the straightforward finite difference ad- 
vection scheme, indicating that this scheme artificially 
enhances the correlation time of the Maxwell stress. For 
the standard finite difference advection scheme of the 
Pencil Code we furthermore observe a tendency for tur- 
bulence to be weaker near the edge of the box and for 
the largest radial density mode to have its minimum near 
the center of the box (see Fig. [16] of Appendix ICl) for the 
two la rgest boxes, confirming the results of Johnso n et al.l 
([2008'). However, with the two improved shearing box 
algorithms Systematic Radial Displacement (SRD) and 
Shear Advection by Fourier Advection (SAFi), any con- 
sistent preference for stronger turbulence near the box 
center is no longer present. This was our primary reason 
for presenting the SAFI versions of L and H in the main 
text. 

7.2. Convergence Tests 

To quantify the effect of dissipation parameters and 
grid resolution on the amplitude of zonal flows and on 
the inverse cascade of magnetic energy we run two sets of 
simulations: (a) runs MJ3z at three different resolutions 
with decreasing diffusivity coefficients and (b) runs MS 
at two different resolutions, the higher resolution simu- 
lations run both with a fixed and with a decreased diffu- 
sivity. The results are given in Table [2] and Table [3l 

For the imposed field simulations M_Bz there is an 
overall increase in the turbulent activity of 50% when 
going from the lowest (M_Bz_r0.5) to the highest reso- 
lution (M_Bz_r2), while the zonal flow and the density 
bumps decrease by around 20% in amplitude. In this 
series of runs we have decreased the hyperdiffusivity co- 
eflScients proportional to Sx^ with increasing resolution. 
The increase in turbulent stresses may thus be due to 
the ever larger inertial range. At the same time these 
added scales will act as an increased diffusivity on the 



large scale pressure bumps, so it is not surprising that 
the zonal flows fall slightly in amplitude. 

The MS series of runs are non-stratified and have 
a computational domain of ~ Ly ~ 2.64_ff and 
Lz — 1.32H. For the low resolution simulation (MS) 
we have forced a time-step similar to that of the higher 
resolution run (MS_r2_fix), to have approximately the 
same amount of numerical dissipation, and we apply the 
Systematic Radial Displacement scheme described in Ap- 
pendix [O to ensure that the numerical dissipation does 
not have significant dependence on x. Wc find that 
zonal flows decrease in amplitude when increasing the 
resolution while keeping the mesh hyper-Reynolds num- 
ber constant (compare run MS to run MS_r2 in Table 
[3l). The mesh hyper- Reynolds number is here defined as 
Rca = max{u)Sx^ /v^. However, the overall turbulent ac- 
tivity also falls significantly as the resolution is increased, 
so the reduction in the large scale amplitudes may sim- 
ply be due to the decrease in available power at the small 
scales. On the other hand if dissipation parameters are 
kept constant as the resolution is increased (as in run 
MS_r2_fix), then there is little or no effect on turbulent 
activity and on zonal flow/pressure bump statistics. It 
is encouraging that extra spatial resolution has so little 
effect on the turbulent state when the dissipation, and 
thus spectral range, is held fixed. Ultimately even higher 
resolution simulations will be needed to fully discern the 
effect of scale separation and Reynolds numbers on the 
zonal flow amplitude. 

7.3. Laplacian Dissipation 

The simulations we have presented so far all used 
hyperdiffusivity operators to dissipate energy near the 
smallest scales of the simulation. This is a "trick" that 
allows us to extend the inertial range of the turbulence, 
at the cost of a steeper transition into the dissipative 
subrange. To check that the use of hyperdiffusivity does 
not create any spurious flows, we have run three simula- 
tions using regular Laplacian operators for viscosity and 
resistivity (runs SS_r4_nul, SS_r8_nul, and MS_r4_iml). 
In order to ensure that the box can sustain turbulence, 
we set the magnetic Prandtl number (defined as ratio of 
viscosity coefficient vi to resistivity coeffi cient 7?i) to ap- 
proxi mately Pm = /r/i = 4 (following iFromang et al.l 
120071 who found that turbulence in zero net flux boxes 
without stratification is generally only sustained when 
Pm > 1). Still the Pencil Code only showed sustained 
turbulence at a resolution of at least 100 grid points per 
scale height, as lower resolution would either show de- 
caying turbulence at low magnetic Reynolds numbers, 
or crash due to lack of dissipation at high magnetic 
Reynolds numbers. 

For moderate and high resolutions, the results in Table 
[2] and Table [3] show that the Laplacian simulations ex- 
hibited fairly strong turbulence with just as strong zonal 
flows as the hyperdiffusivity runs. This yields much sup- 
port to the fact that the zonal flows are real, and not a 
numerical artifact, because of the more physical nature 
of the Laplacian operators and the fact that zonal flows 
operate also at 200 grid points per scale height in res- 
olution, a factor eight times higher resolution than the 
hyperdiffusivity runs S-H. However there is a significant 
increase in turbulent activity when going from run SS_r4 
to run SS_r8. The fact that the MRI has not yet con- 
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TABLE 4 

Zonal Flow Properties for Various Shearing Box Schemes 



Run 


{ — Bx By) 




\Uy\{k0,0,0) 


|p|(fco,0,0) 




|S^|(fco,0,0) 






Cs B 

y 


^corr 


S (FDA) 


4.2 X 10" 


3 


1.9 X 10" 


2 


8.1 X 10" 


3 


4.5 X 10 


-4 


-1.6 


-2.4 


-0.9 


8.3 


S (RRC) 


2.3 X 10" 




1.8 X 10" 


2 


8.0 X 10" 




2.7 X 10" 


-4 


-1.5 


—2.4 


-0.7 


14.9 


S (SRD) 


3.1 X 10" 




1.8 X 10" 


2 


7.9 X 10" 




3.8 X 10" 


-4 


-1.5 


-2.3 


-0.8 


8.9 


b (SAl^l) 


2.8 X 10" 


3 


1.8 X 10" 


2 


8.1 X 10" 


3 


3.7 X 10 


-4 


— 1.5 


—2.5 


—0.8 


12.3 


M (FDA) 


1.2 X 10" 


■2 


2.4 X 10" 


2 


1.9 X 10" 


2 


1.2 X 10 


-3 


-1.1 


-1.8 


-0.7 


8.5 


M (RRC) 


8.6 X 10" 


3 


2.6 X 10" 


2 


2.4 X 10" 


2 


9.7 X 10 


~4 


-1.3 


-2.3 


-1.0 


11.6 


M (SRD) 


8.8 X 10" 


3 


1.7 X 10" 


2 


1.7 X 10" 


2 


8.7 X 10 


-4 


-0.5 


-1.8 


-0.8 


8.3 


M (SAFI) 


9.6 X 10" 


3 


1.9 X 10" 


2 


1.6 X 10" 


2 


9.2 X 10 


-4 


-0.8 


-1.8 


-0.7 


8.4 


L (FDA) 


1.1 X 10" 


2 


2.5 X 10" 


2 


5.0 X 10" 


•1 


1.3 X 10 


-3 


-1.1 


-2.3 


-1.0 


18.6 


L (RRC) 


1.3 X 10" 


2 


1.8 X 10" 


2 


3.6 X 10" 


2 


1.1 X 10- 


-3 


-0.5 


-1.6 


-0.5 


20.1 


L (SRD) 


1.1 X 10" 


2 


2.0 X 10" 


2 


4.0 X 10" 


2 


1.2 X 10- 


-3 


-0.6 


-1.9 


-0.9 


20.9 


L (SAFI) 


1.1 X 10" 


2 


1.7 X 10" 


2 


3.3 X 10" 


2 


1.3 X 10- 


-3 


-1.0 


-2.3 


-1.4 


19.6 


H (FDA) 


9.5 X 10" 


3 


1.6 X 10" 


2 


5.8 X 10" 


•2 


1.1 X 10 


-3 


-0.4 


-1.4 


-1.1 


57.1 


H (SAFI) 


1.0 X 10" 


2 


1.5 X 10" 


2 


5.2 X 10" 


2 


1.1 X 10 


-3 


-0.2 


-1.1 


-0.9 


50.5 



Note. — Col. (1): Name of run. Col. (2): Maxwell stress, normalized by mean pressure. 
Col. (3)- (5): Fourier amplitude of azimuthal velocity, density normalized by mean density in the 
box, and pressure- normalized Maxwell stress, respectively, at largest radial mode. Col. (6)-(8): 
Power law index of azimuthal velocity, density, and Maxwell stress, respectively, calculated 
from two largest radial modes. Col. (9): Correlation time, in orbits T = 2-kO~^, of largest 
radial density mode. The various shearing box algorithms - FDA, RRC, SRD, and SAFI - are 
defined in Appendix [Cl 



verged even at such high resolution is a slight worry. It 
may indicate either that numerical dissipation still influ- 
ences the saturation level of the turbulence or that the 
introduction of smaller scales changes the dynamics of 
the system. One would hope that the MRI will converge 
at the next resolution level, i.e. 400 grid points per scale 
height, but this huge computational effort is beyond the 
scope of the current investigation. The ratio of the large 
scale Maxwell stress variation to the mean value never- 
theless stays approximately the same for the two runs 
SS_r4 and SS_r8, indicating some convergence in the in- 
verse cascade of magnetic energy. 

8. SUMMARY AND DISCUSSION 

We conduct numerical simulations of turbulence driven 
by the magnctorotational instability (MRI) in large 
shearing boxes. We here emphasize four main results. 
(1) The turbulent energy and stresses more than dou- 
ble when the box size (radial and azimuthal) increases 
from 1.32 to 2.64 scale heights. In yet larger boxes the 
turbulent fluctuations and stresses remain relatively con- 
stant. (2) Magnetic energy exceeds kinetic energy by a 
factor two-three on scales above approximately two scale 
heights, when the box size is sufficient for the inverse 
case to proceed to these large scales. (3) Axisymmetric 
surface density fluctuations grow to fill the box, and per- 
sist for many orbits, with an increasing lifetime in bigger 
boxes. These pressure bumps are in gcostrophic balance 
with sub/super-Keplerian zonal ffows. (4) We put for- 
ward a simple model for the launching of zonal flows by 
large scale variations in the Maxwell stress. This ran- 
dom walk model explains both the saturated amplitude 
as well as the correlation time of pressure bumps and 
their associated zonal flows. 

Large scale pressure and zonal flow structures grow to 
fill the simulation domain for all box sizes. The largest 
box size considered, more than ten scale heights in ra- 
dial and azimuthal extent, nevertheless indicates an in- 



creasingly flat dependence of zonal flow power on the 
radial wave number. Even larger simulation domains 
will be needed to verify whether the natural scale of the 
large scale structures is really around ten scale heights. 
Global effects - including curvature and radial variation 
in background quantities like density and sound speed 
- may also play a vital role for the termina tion of the 
invers e cascade. The global simulations by iLvra et all 
(1200^ show evidence of long lived radial pressure bumps 
that have not cascaded to the largest radial scale of the 
disk (see their Fig. 18). It is intriguin g to note that the 
length scale of the pressure bumps in iLvra et all (|2008l ) 
is around ten scale heights. 

We have experimented with several improvements of 
the shearing box algorithm in order to quantify the de- 
pendence of the results on numerical diffusivity associ- 
ated with the Kcplcrian shear ffow. There is a ten- 
dency for the density depression to occur at the cen- 
ter of the box for the two largest box sizes (approxi- 
mately five and ten scale heights in radial extent, re- 
spectively) if shear advection is done by the usual finite 
difference scheme of the Pencil Code, co nfirming a simi- 
lar eff ect found bv I Johnson et all (|200l for a ZEUS-hke 
code (|Stone fc Normanlll992l )~ However, when shear ad- 
vection is integrated by high order interpolation instead, 
we see no preference for any special location in the box. 
The same is true if the state of the gas is displaced sys- 
tematically by one grid point in the radial direction at 
the beginning of each time-step. 

Ultimately the emergence of zonal fiows and pres- 
sure bumps is controlled by the large scale variation 
of the magnetic field. The cascade of magnetic energy 
from small to large scales only occurs for the x- and y- 
components of the field. The vertical field, which is of 
importance for feeding the vertical field MRI, shows very 
little large scale behavior. The specific mechanism that 
leads to the inverse cascade of the in-plane field in our 
simulations is not known. In the absence of any net ki- 
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netic helicity separation in the radial direction, the large 
scale field may ari se as an effect of a negative turbu- 
lent diffusivity (e.g. iZheligovskv et al.ll200ll: |Ur pin"2002f) 
or a shear current effect (jRogachevskii fc Klec orin 2003; 
lYousef et al.][200l . 

Both hyperdiffusivity simulations and direct numerical 
simulations (DNS), with Laplacian viscosity and resis- 
tivity operators, show prodigious zonal flows. A fruit- 
ful future path could be to focus on even higher res- 
olution direct simulations of resistive MRI turbulence. 
This could probe the dependence of the inverse cascade 
of magnetic energy and the launching of zonal flows on 
the magnetic Reynolds number and the Prandtl num- 
ber (|Lesur fc Longaretti|[2007HFromang et al.ll2007D . be- 
lieved to be m uch smaller than unity in most parts of 
accretion disks (jBalbus fc Henril 120081 ). 

Our numerical experiments show that the magnetic 
field takes particular advantage of larger spatial scales 



by participating in an inverse cascade. This in turn pro- 
duces long-lived density fluctuations which can have sig- 
niflcant influence on the formation and evolution of plan- 
etary systems. Big shearing boxes thus provide an excel- 
lent tool for studying turbulent accretion disks and the 
physical processes that occur therein. 



We are grateful to Takayoshi Sano, Ethan Vishniac, 
Axel Brandenburg, Eric Blackman, Yuri Levin and 
Wladimir Lyra for inspiring discussions on large scale 
magnetic fields in accretion disks. We would like to thank 
the anonymous referee for an insightful referee report. 
Bryan Johnson is thanked for commenting an early ver- 
sion of the manuscript. H. K. has been supported in part 
by the Deutsche Forschungsgemeinschaft DFG through 
grant DFG Forschergruppc 759 "The Formation of Plan- 
ets: The Critical First Growth Phase" . 



REFERENCES 



Arlt, R., & Riidigcr, G. 2001, A&A, 374, 1035 
Armitagc, P. J. 1998, ApJ, 501, L189+ 
Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 21 
Balbus, S. A., & Hawley, J. F. 1992, ApJ, 400, 610 
Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physies, 
70, 1 

Balbus, S. A., & Hawley, J. F. 2006, ApJ, 652, 1020 

Balbus, S. A., & Henri, P. 2008, ApJ, 674, 408 

Bartello, P., Metais, O., & Lesieur, M. 1994, Journal of Fluid 

Mechanics, 273, 1 
Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 
Borue, V. 1993, Physical Review Letters, 71, 3967 
Brandenburg, A., Nordlund, A., Stein, R.F., & Torkelsson, U. 

1995, ApJ, 446, 741 
Brandenburg, A. 2003, in Advances in nonlinear dynamos (The 

Fluid Mechanics of Astrophysics and Geophysics, Vol. 9), ed. 

A. Ferriz-Mas & M. Nuricz (Taylor & Francis, London and New 

York), 269-344 

Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1 
Brauer, F., DuUemond, C. P., & Henning, T. 2008, A&A, 480, 859 
Brauer, F., Henning, T., & DuUemond, C. P. 2008, A&A, 487, LI 
Busse, F. H. 1976, Icarus, 29, 255 

Cambon, C., Mansour, N. N., & Godeferd, F. S. 1997, Journal of 

Fluid Mechanics, 337, 303 
Diamond, P. H., Itoh, S.-L, Itoh, K., & Hahm, T. S. 2005, Plasma 

Physics and Controlled Fusion, 47, 35 
Dominik, C, Blum, J., Cuzzi, J. N., & Wurm, G. 2007, in 

Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & 

K. Keil, 783-800 
DuUemond, C. P., & Dominik, C. 2005, A&A, 434, 971 
Foglizzo, T., & Tagger, M. 1995, A&A, 301, 293 
Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in 

Astrophysics: Third Edition (Cambridge University Press, 

2002) 

Fromang, S., & Nelson, R. P. 2005, MNRAS, 364, L81 
Fromang, S., & Nelson, R. P. 2006, A&A, 457, 343 
Fromang, S., & Papaloizou, J. 2007, A&A, 476, 1113 
Fromang, S., Papaloizou, J., Lesur, G., & Heinemann, T. 2007, 

A&A, 476, 1123 
Gammie, C. F. 2001, ApJ, 553, 174 
Goodman, J. 2003, MNRAS, 339, 937 
Haghighipour, N., & Boss, A. P. 2003, ApJ, 583, 996 
Haugen, N. E., & Brandenburg, A. 2004, Phys. Rev. E, 70, 036408 
Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 

742 

Heimpel, M., & Aurnou, J. 2007, Icarus, 187, 540 
Hirose, S., KroUk, J. H., & Stone, J. M. 2006, ApJ, 640, 901 
Howard, R., & Labonte, B. J. 1980, ApJ, 239, L33 
Ida, S., Guillot, T., & MorbideUi, A. 2008, ApJ, 686, 1292 
Johansen, A., & Klahr, H. 2005, ApJ, 634, 1353 
Johansen, A., Klahr, H., Henning, Th. 2006, ApJ, 636, 1121 
Johansen, A., Oishi, J. S., Low, M.-M. M., Klahr, H., Henning, 
T., & Youdin, A. 2007, Nature, 448, 1022 



Johansen, A., Brauer, F., DuUemond, C, Klahr, H., & Henning, 

T. 2008, A&A, 486, 597 
Johnson, E. T., Goodman, J., & Menou, K. 2006, ApJ, 647, 1413 
Johnson, B. M. 2007, ApJ, 660, 1375 

Johnson, B. M., Guan, X., & Gammie, C. F. 2008, ApJS, 177, 373 
Kato, M. T., Nakamura, K., Tandokoro, R., Fujimoto, M., & Ida, 

S. 2008, ArXiv e-prints, 0810.3466 
Klahr, H., & Lin, D. N. C. 2001, ApJ, 554, 1095 
Kraichnan, R. H. 1967, Phys. Fluids, 10, 1417 

Laughlin, G., Steinacker, A., & Adams, F. C. 2004, ApJ, 608, 489 
Lesur, G., & Longaretti, P.-Y. 2007, MNRAS, 378, 1471 
Lesur, G., & Ogilvie, G. I. 2008, A&A, 488, 451 
Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2008, A&A, 
479, 883 

Masset, F. 2000, A&AS, 141, 165 

Miller, K. A., & Stone, J. M. 2000, ApJ, 534, 398 

Moffatt, H. K. 1978, Magnetic field generation in electrically 

conducting fluids (Cambridge University Press, 1978) 
Nelson, R. P., & Papaloizou, J. C. B. 2004, MNRAS, 350, 849 
Oishi, J. S., Mac Low, M.-M., & Menou, K. 2007, ApJ, 670, 805 
Parker, E. N. 1955, ApJ, 122, 293 
Pringle, J. E. 1981, ARA&A, 19, 137 

Rogachevskii, I., & Klecorin, N. 2003, Phys. Rev. E, 68, 036301 
Safronov, V. S. 1969, Evoliutsiia doplanetnogo oblaka. (English 

transl.: Evolution of the Protoplanetary Cloud and Formation 

of Earth and the Planets, NASA Tech. Transl. F-677, 

Jerusalem: Israel Sci. Transl. 1972) 
Schiissler, M. 1981, A&A, 94, L17 
Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 
Sheehan, D. P., Davis, S. S., Cuzzi, J. N., & Estberg, G. N. 1999, 

Icarus, 142, 238 
Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 791 
Sun, Z.-P., Schubert, G., & Glatzmaier, G. A. 1993, Science, 260, 

661 

Terquem, C, & Papaloizou, J. C. B. 1996, MNRAS, 279, 767 

Urpin, V. 2002, Phys. Rev. E, 65, 026301 

Vishniac, E. T., & Brandenburg, A. 1997, ApJ, 475, 263 

Vishniac, E. T., submitted to ApJ 

Weidenschilling, S. J. 1977, MNRAS, 180, 57 

Weidenschilling, S. J. 1997, Icarus, 127, 290 

Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211 
Yoshimura, H. 1981, ApJ, 247, 1102 
Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494 
Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 
Youscf, T. A., Heinemann, T., Schekochihin, A. A., Klecorin, N., 
Rogachevskii, I., Iskakov, A. B., Cowley, S. C, & McWilliams, 
J. C. 2008, Physical Review Letters, 100, 184501 
ZheUgovsky, V. A., Podvigina, O. M., & Frisch, U. 2001, 
Geophysical and Astrophysical Fluid Dynamics, 95, 227 



Zonal Flows in Magnctorotational Turbulence 



19 




Fig. 12. — Normalized kinetic energy of a time-dependent non-axisymmetric shear wave, for different numerical resolutions and the 
analytical solution. The numerical solution follows the analytical solution well down to around 6-8 grid points per radial wavelength, after 
which the explicit viscosity is so strong that the shear wave dies out. There is no evidence of non-linear aliasin g as the wave approaches 
the grid scale, again due to damping by the explicit viscosity terms at the small scales. Compare with Fig. 3 of lBalbus &: Hawlevl I I2006I) 
which shows the results with the Athena code. 

APPENDIX 

A. TEST OF SHEARING BOX IMPLEMENTATION 

In this Appendix we test the shearing box implementation of the Pencil Code against known analytical solutions to 
the evolution of shearing waves. In i ^A.ll we compare the results of the Pencil Code to the analytical solution fou nd by 
iBalbus fc Hawlevl (|2006f ) for hydrodynamical shearing waves. We proceed to test the magnetic field module in i]A.2i 
where we compare a numerical solution of the linearized shearing box MHD equations to the results obtained with the 
Pencil Code. We find excellent agreement in both cases. 

A.l. Time- Dependent Incompressible Shear Wave 

IBalbus fc Hawlevl (|2006f ) derived a time-dependent, non-axisymmetric and incompressible solution for hydrodynam- 
ical shear waves. Solutions are of the form of a single mode f{x,y,z,t) = f(t)cxjp[i{kx{t)x + kyy + k^z)] for the 
dynamical variables u^, Uy and u^. The Keplerian shear causes the radial wavenumber to be time-dependent and have 
the form k^^it) = k^{0) + l?,/2)Qtky. 

We have performed numerical simulations of the model set up described in §5 of lBalbus fc Hawlevl (|2006D . We 
use the simulation parameters q = 3/2, Q = 0.001, p = 1, Cg = 4.08^?, 7 = 5/3 and P = 10"^. The box size is 
Lx = Ly = Lz = 1. The initial condition is a non-axisymmetric wave in the x- velocity, 

Ux{x, y, t = 0) = Ar sin[27r(y + Az)] , (Al) 

with the amplitude Ar ~ 10^^. 

In Fig. [T2|we show the numerical solution for a range of resolutions and compare with the analytical solution. There 
is excellent agreement between the Pencil Code integration and the analytical solution, down to 6-8 grid points per 
wavelength where the viscous damping is so stro ng that the wave d i es out . There is no trace of non- linear aliasing as 
the wave is sheared below the spatial resolution. IBalbus fc Hawlevl (|2006D show a similar plot for the Athena code in 
their Fig. 3. Actually the Pencil Code is slightly less dissipative than the Athena code, but this of course depends on 
our choice of dissipation parameters. For the results in Fig.[T2|we set the hyperviscosity and hyperdiffusion similar to 
their values in the MRI simulations (i.e. similar viscous dissipation time at the grid scale). 

A. 2. Magnetized Shear Wave in the Linear Regime 

To test the evolution of the magnetic vector potential in the shearing box we set up a magnetized shearing wave 
test problem. Magnetized shear waves are an interesting analytical and numerical topic, but in this app endix we shall 
consider such waves solely to test the numerical solver of the Pencil Code. We refer to iJohnsonI (|2007D for details on 
the physical behavior of magnetized shear waves. 

Linearizing the equation of motion, continuity equation and induction equation as / = /o + /' for all components 
and considering the time evolution of individual Fourier modes of the form 



f'{x, y, z, t) = fit) exp[i(fc2;(t)x + kyy + k^^z)] , 



(A2) 
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Fig. 13. — Kinetic and magnetic energy of a magnetized, non-axisymmctric shear wave as a function of radial wavenumber (left plots). 
The relative error between analytical and numerical solution is shown in the right plots (here zero is defined as no error, while an error 
measure of unity corresponds to infinite error). There is excellent agreement between the numerical solution of the Pencil Code and the 
analytical solution, down to approximately 6 grid points per wavelength. 



with time-dependent radial wavenumber 



yields 



— - Uye^-- u^ey + — 

— = -poifc • u , 

3 

— =i(Bo ■k)u- -QB^ey - iBak-ii. 



-ik{Bo ■ B)+i{Bo ■ k)B 



— ife/5, 



(A3) 

(A4) 
(A5) 
(A6) 



Setting B = {0,0, Bq), k — (0,0, fc^) ii = (ux,Uv,0) and B = {Bx,By,0) would yield evolution according to the 
linear magnetorotational instability (|Balbus fc HawlevI |1991|) . But since we are interested in the potential effect 
of shear advection and shear-periodic boundary conditions we instead seed a compressive, non-axisymmetric mode 
with k = [-2,1,4], ii = [0.001,0,0], po = 1, J5o = [0,1,0], B = [0,0,0] in a box of size (27r)3. Such a leading 
wave perturbation of an azimuthal field line is subject to a transient non-axisy mmctric magnetorotational instability 
(IBalbus fc Hawlevl[l99l fFoglizzo fc Tagged [l995l : iTerauem fc PapaloizoulllOOlh . In Fig. [H we compare the evolution 
of magnetic and kinetic energy obtained by the Pencil Code to the evolution obtained by integrating the linearized 
equation system in time. There is excellent agreement between the two down to approximately 6 grid points per 
wavelength. 
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B. KEPLERIAN SHEAR ADVECTION AS INTERPOLATION 

In this Appendix we describe the numerical implementation of interpolated Keplerian advection in the Pencil Code. 
The spl itting of the ad vection step into advection by Keplerian she ar and a dvect ion by the fluctuating velocity goes 
back to lGammi^ ()2001[ ) for the shearing box approximation and to iMasseO (|2000( ) for global, cylindrical coordinates. 
The interpolation scheme was recentl y generalized to inclu de magnetic fields, in the shearing box approximation, and 
implemented in a ZEUS-like code by Ijohnson "eFall ((2OOI . 

All shearing box dynamical variables / have a Keplerian advection term of the form 

dt y dy ^ ^ 

in their time evolution equation. Here u^y"^ = —{Z/2)flx is the linearized Keplerian velocity. The time-step associated 

with the Keplerian advection is 5t = cst^x / v[i&y^{\u"y\) , where est is a factor of less than unity (in the Pencil Code 
stability requires approximately est < 0.4). The Keplerian advection time-step dominates over the time-step from 

sound waves when max(|My°''|) > Cg, i.e. when Lx/H > 2/(3/2) w 1.33 for Keplerian rotation. 
The time-step of the Keplerian advection is irrelevant if that term is not treated by finite-differencing, but rather as 

a shift of all variables in physical space by the amount u^^St. Without deriving it formally we can estimate that the 
stability criterion with this method is that neighboring points in the ai-direction must not move more than a fraction 
^ of a grid point apart, 

{3/2)nSxSt<^Sy. (B2) 

This is a much less tight time-step constraint than for the finite difference approach, since there is no longer any 
reference to the radial coordinate x in the time-step calculation. The interpolation scheme also eliminates the numerical 
diffusivity associated with the Keplerian advection (see Appendix [Ujl . 
We implement the interpolated Keplerian advection as follows: 

1. All variables q{x, y, z) are Fourier transformed in the j/-direction to yield q{x, ky, z). 

2. Each Fourier amplitude is multiplied by the complex factor eKp[ikyU^\x)dt\ to shift by the amount u^y\x)5t in 
real space. 

3. Inverse Fourier transform to real space. 

Because of the third order Runge-Kutta time-integration scheme of the Pencil Code, the shift must be done for the 
inherited right-hand-sides q at the beginning of the second and third sub-time-steps as well. The time-step entering 
in point 2 above, now dominated by sound waves, is calculated at the first sub-time-step of the Runge Kutta scheme 
and consequen tly used for all three substeps, including the Keplerian advection. 

Note that in [Massed (j200(3( ) the state is always shifted by an integer number of grid points, treating the remaining 
Keplerian advection speed together with the perturbed speed. In this paper, on the other hand, the entire Keplerian 
advection is done by interpolation in Fourier space. We denote this scheme Shear Advection by Fourier Interpolation 
(SAFI). 

C. REMOVING POSITION-DEPENDENT NUMERICAL DIFFUSIVITY 

In this Appendix we quantify the numerical diffusivity of the Pencil Code advection scheme and, very importantly 
for simulations of turbulence in large shearing boxes, its dependence on the background advection speed. We show 
how to eliminate the spatial dependence of the numerical diffusivity by various improvements to the shearing box 
algorithm of the Pencil Code. 

C.l. Measuring the Numerical Diffusivity of the Pencil Code 

Although numerical diffusion is not necessarily representable as a well-behaved differential operator, we shall fit the 
time evolution of a passively advected physical variable / according to the expression 

|/|(t) = 1/1(0) exp(-i?„fc2"t). (CI) 

The hat denotes here that / is set as a pure Fourier mode at scale k, and J5„ is a diffusivity cocfRcicnt of order n. 
Regular diffusion would e.g. have n = 1, while hyperdiffusion has higher values of n. 

We measure the diffusivity of the Pencil Code's advection scheme by advecting a wave of arbitrary initial amplitude 
for the time Ai = 100 with the constant advection speed uq. We then fit the amplitude reduction parameter, 

Q^J^im). (C2) 

I, 1/1(0) / 

according to the power law 

Qcxfc2»u™(Ji)'. (C3) 
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Fig. 14. — The dependence of the amphtude reduction parameter Q = Dnk^"At on wavenumber k (first plot), advection speed uq (middle 
plot) and time-step St (last plot). All dependencies are fitted excellently by a power-law, indicating an numerical diffusivity that behaves 
like second order hypcrdiffusivity. 

Here St is the time-step of the temporal integration. If the wavenuniber-dependence can indeed be fitted with a single 
n, then the diffusion coefficient is fitted according to £)„ oc it™((5t)'. We show in Fig. [T3]the dependence of Q on k, 
uq and St. There is an extremely good power-law fit for each parameter, with n = 2, m = 4 and I = 3. Thus the 
numerical diffusivity behaves like a second order hyperdiffusion. 



with diffusion coefficient 



D2 ~ OMlu'^idt)^ 



(C4) 



(C5) 

One may in fact show from stability analysis that D2 ~ {1 / 24:)uQ{St)^ (W. Dobler, private communication). The 
strong dependence on the time-step St is not surprising for the 3rd order Runge-Kutta scheme of the Pencil Code. One 
consequence is that if St is determined from the global maximum speed, and uq varies across the simulation domain, 
then the diffusion coefficient goes as the background speed to the fourth power, which can lead to huge differences in 
the numerical diffusivity across the box. 

C.2. Improved Shearing Box Algorithms 

We proceed to test the effect of the varying numerical diffusivity in a linear shear fiow. We have set up a simple 
advection test where a wave of azimuthal wavelength 8 grid points is advected by a linear shear fiow with velocity 
difference (3/2)7r « 4.71 from box edge to box center (the radial extent of the box is = 2tt). We follow the amplitude 
of the wave as a function of the distance from the center of the box. This is plotted in Fig. [151 The first panel shows 
that for the usual advection scheme the error grows much quicker near the box edges, while the box center is left 
almost untouched by numerical diffusion^. 

There are several ways to reduce, or even remove, this space-dependent numerical diffusivity. Here we discuss three 
ways: 

1. Picking the center of revolution randomly at each time-step, thus making sure that even the central points are 
advected by the shear. 

2. Displacing the entire data cube by one grid point in the radial direction in the beginning of each time-step. 

3. Treating the shear advection term separately by high order interpolation. 

C.2.1. Random Revolution Center (RRC) 
The linear shear flow can be generalized as 



(C6) 



where xq is the radial coordinate of the revolution center, normally set to zero. We can nevertheless choose an arbitrary 
xo G [—L,j;/2, Lx/2] at each time-step and still solve exactly the same equation system. This way we can reduce the 



The central grid points i = 3,4 lie at both sides of x = 0, thus even these grid points experience some shear advection. 
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Fig. 15. — The amplitude of a passively advectcd shear wave as a function of time t at four radial positions in the box (different colors). 
The amplitude falls due to the numerical diffusivity associated with the Keplerian advection term. The simple finite difference scheme gives 
much higher error near the edges of the box (first plot), as the numerical viscosity scales with the local advection speed to the fourth power. 
Centering the shear around a random point in the box increases the error in the center (second plot), but the error maintains a significant 
dependence on the radial coordinate. By displacing the entire data cube by one grid point in the radial direction in the beginning of each 
time-step or by interpolating the shear advection term in Fourier space one can remove the radial dependence of the amplitude error (third 
and fourth plots). 

special state that a; = normally has because that point is not advected. The second plot of Fig. [T5l shows the radially 
dependent amplitude error using the Random Revolution Center scheme. While the error near the center is greatly 
increased compared to the usual advection scheme, there is still a significant difference in the amplitude between the 
center and the edges of the box. 

The time-step constraint of the RRC method is slightly tighter than for the usual advection scheme, since the 
maximum advection speed is around a factor two times higher. The higher average advection speed leads to some 
increase in the overall numerical diffusivity, although the lower time-step partially counteracts this increase. 

C.2.2. Systematic Radial Displacement (SRD) 

A better way to reduce the radially dependent numerical diffusivity is to distribute the error over the entire box. We 
do this by shifting the data cube by one grid point in the radial direction at the beginning of each (main) time-step. 
In the time interval At = Nx5t a physical structure will have experienced numerical diffusivity at all radial locations 
in the box, and since this is much shorter than the relevant physical time-scale 27ri7~^, the total numerical diffusivity 
should be approximately the same for all flow features. In the third plot of Fig. [15] we show how all grid points have 
the same amplitude error, even after a long integration time. The time-step may furthermore be lowered artificially 
in order to reduce the time spent at any given radial location. This method is preferred over the Random Revolution 
Center scheme, but at the cost of a slightly higher numerical diffusivity as the shearing radial boundary conditions 
must be set twice each time-step. Also postprocessing and on-the-fly diagnostic output is complicated by the rapid 
radial displacement of the grid. 

C.2.3. Shear Advection by Fourier Interpolation (SAFI) 

A third possibility is to reduce amplitude error in the advection scheme through time integration of the Keplerian 
advection term by interpolation. The Shear Advection by Fourier Interpolation scheme is described in detail in 



24 



Johansen, Youdin, & Klahr 



FDA 



RRC 



SRD 



SAFI 



s 


s 


s 


S ■ 




ll 






M 


M 


M 


M ■ 






A^^^ ^J^T-^^^j-H^ 






L 


L 

^ ^»*j^yi]. 


L 

LhT \jln.^^ m 


L ■ 








H ■ 



c 


5 

ej 
c 

05 
i3 

^3 



;.o 2- 

.5 

n 

1.0 

;.o2- 

.5 2 

.ol 

1.5.2 
Q 

1.0 



-4.0 -2.0 0.0 2.0 4.0 
x/H 



-4.0 -2.0 0.0 2.0 4.0 
x/H 



Fig. 16. — The Maxwell stress, divided by the mean Maxwell stress, as a function of the radial distance from the center of the box (black 
curve). Overplotted is the distribution function for the location of the minimum of the largest radial density mode (gray histograms). Rows 
represent various box sizes, while columns show the results for Finite Difference Advection (FDA), Random Revolution Center (RRC), 
Systematic Radial Displacement (SRD), and Shear Advection by Fourier Interpolation (SAFI). The typical fluctuation in the Maxwell 
stress is in all cases around 10% (see Table [Sj. In the smallest boxes these fluctuations have almost canceled out each other, since the 
correlation time of the stress is so short. The two largest boxes, on the other hand, have an excess of a few percent. Here the 100 orbits 
integration has not covered enough coherence time-scales to give a good average. There is a slight tendency for the turbulence to be weak 
near the edges of the box, and for the density minimum to appear near the center of the box, when using the standard finite difference 
advection scheme (FDA, first column). 

Appendix [B] where we focused on making to allowed time-step longer by integrating the Keplerian advection term 
through high order interpolation. This scheme however also reduces the amplitude error of the advection to zero, 
yielding a numerical difFusivity that does not depend on the radial coordinate. The last panel of Fig. [T5] shows that 
the amplitude error is indeed zero with this scheme. 

In fact the SAFI method advects all scales with zero amplitude error and zero dispersion error. The only exception 
is the Nyquist scale, which is damped as the wave approaches a cosine function of y. This process of damping power 
at the Nyquist scale can have a stochastic element to it, as it may take arbitrarily long for the sine wave to be mapped 
as a cosine wave on the grid. 

So in effect there will be an x-dependent amplitude error at the Nyquist scale. This is however unlikely to affect 
the overall diffusivity of the Keplerian advection, since the Nyquist scale is already heavily damped by the explicit 
(hyper, shock, or regular) diffusivity terms. A combination of the SRD and SAFI schemes would in theory remove 
the .T-dcpendcnt damping of the Nyquist scale, but we have not explored this any further since the SRD scheme 
complicates postprocessing significantly. 



C.3. Comparison of Results Obtained with Various Shearing Box Schemes 

To test the relevance of space-dependent numerical dissipation for our results we have run four variations of each 
of the standard box sizes S, M, L (FDA, RRC, SRD, SAFI), and two variations of the largest box H (FDA, SAFI). 
The Maxwell stress, normalized by the mean Maxwell stress, is shown in Fig. [11] as a function of radial distance from 
the center of the box. We have overplotted the distribution function for the location of the minimum of the largest 
radial density mode (see histograms and the right axis). For the standard finite difference advection scheme there 
is a tendency for the turbulence to be weaker near the edges of the box and for the density minimum to occur near 
the center of the box (see first column). However, the preference for turbulence to be strong near the box center is 
not consistently present with any of the improved shearing box algorithms. Thus we conclude that space-dependent 
numerical diffusivity can potentially lead to spurious flow features for large shearing boxes with the Pencil Code, but 
that this problem is alleviated by either displacing the grid systematically in the radial direction at each time-step or 
by integrating the Keplerian advection term by Fourier interpolation. 

The measured statistical properties of zonal flows, pressure bumps and large scale power for various shearing box 
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algorithms arc written in Table 21 For each box size there is a tendency for the large scale power in zonal flows and 
pressure bumps to decrease when going from the standard finite difference advection scheme to the Fourier interpolation 
method. The amplitude of the Maxwell stress, on the other hand, is relatively independent of the chosen algorithm. 
Based on these measurements we found it safer to use the Fourier interpolation scheme for box sizes L and H. 



